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Abstract 

Among various algorithms designed to exploit the specific properties of quan- 
tum computers with respect to classical ones, the quantum adiabatic algorithm 
is a versatile proposition to find the minimal value of an arbitrary cost func- 
tion (ground state energy). Random optimization problems provide a natural 
testbed to compare its efficiency with that of classical algorithms. These prob- 
lems correspond to mean field spin glasses that have been extensively studied 
in the classical case. This paper reviews recent analytical works that extended 
these studies to incorporate the effect of quantum fluctuations, and presents 
also some original results in this direction. 
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1. Introduction 

A central issue in computer science is the classification of the difficulty of 
computational tasks, i.e. the existence or not of algorithms with small require- 
ments (in terms of the time of execution and the necessary memory) that per- 
form a given task [1, 2]. This classification is called computational complexity 
theory. A rough distinction between easy and hard tasks is made by distin- 
guishing algorithms that need to perform a number of elementary operations 
growing either polynomially or exponentially with the size of their input. One of 
the central tasks analyzed in this context concerns combinatorial optimization 
problems [3]: given a cost function defined on N variables, each taking a finite 
number of values, the question is to classify families of cost functions such that 
algorithms can, or cannot, find their global minimum by executing a number 
of operations smaller than some polynomial of N. The current consensus is 
that there exist families of cost functions such that no algorithm can achieve 
this goal (this is the famous P^NP conjecture). One example of such difficult 
problems is the graph g-coloring (for q > 3): given a graph, i.e. a collection of 
N vertices and M edges linking some of the pairs of vertices, the cost function 
associates to each of its g-coloring (one out of q colors is chosen for each vertex) 
the number of monochromatic edges, linking two vertices of the same color. 

It was understood above that the term "elementary operation" meant some 
simple process like adding two numbers, or other arithmetic tasks, which can 
all be reduced to logical operations on boolean variables. In that context the 
basic elements of a computer are bits that behave "classically" , i.e. they are in 
a well defined state or 1, that is altered deterministically by logical operations 
involving one or a few of them. But what happens if the basic elements of a 
computer behave "quantumly" , i.e. if instead of bits one deals with "qubits" 
that can not only be in the states or 1 but in any linear combination of 
the two? Will such a quantum computer be able to solve efficiently some of 
the tasks on which classical computers get stuck? These questions were first 
raised in the eighties by Feynman [4] and Dcutsch [5] and opened the way to a 
new branch of science at the interface between computer science and physics, 
known today under the names of quantum computing and quantum information 
theory [6, 7, 8]. 

Several specific quantum algorithms have been discovered since then [5, 9, 
10, 11, 12], providing "quantum speedup" with respect to their fastest classical 
counterparts. Most of them concern arithmetic problems, notably Shor's algo- 
rithm for factoring integers [11], yet none of them solves efficiently a representant 
of the classically hardest problems (the so-called NP-complete ones). A quantum 
analog of the computational complexity theory has been developed [13, 14, 15], 
with the introduction of complexity classes of easy and hard problems, the no- 
tion of difficulty being now with respect to the number of required operations 
on a quantum, instead of classical, computer. 
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A generic strategy to solve optimization problems with a quantum computer, 
called quantum annealing or quantum adiabatic algorithm [16, 17, 18, 19, 20], 
proceeds in the following way (see [21, 22, 23, 24] for reviews of this procedure). 
The evolution of the state of a quantum computer obeys Schrodingcr equation, 
with a time-evolving Hamiltonian H{t) that is controlled by the programmer. 
If in an initialization step the system is prepared in the ground state of a simple 
Hamiltonian iJj, and if the time evolution of the Hamiltonian is slow enough, the 
adiabatic theorem [25] ensures that the system remains, with high probability, 
in the instantaneous ground state at all subsequent times. This property can be 
exploited by driving the Hamiltonian towards one that corresponds to the cost 
function of the optimization problem to be solved, let us call it Hf. Indeed its 
ground state, which is the final state of the system according to the adiabatic 
theorem, provides precisely the answer to the combinatorial optimization prob- 
lem. The crucial question of the efficiency of such an algorithm reduces thus to 
a criterion for the validity of the adiabatic approximation. Roughly speaking 
the adiabatic theorem states that the evolution time of the Hamiltonian (hence 
the running time of the algorithm) has to be larger than the inverse square of 
the minimal energy gap between the ground state and the first excited state 
encountered during the process. Instances of optimization problems Hf such 
that this gap is exponentially small in the size of the problem (and thus require 
an exponentially large time to be solved adiabatically) were exhibited early 
on [26, 27]. It was also realized that some choices of the initial Hamiltonian iJj 
led ineluctably to exponentially small gaps [28, 29]. 

One can however wonder if, for "reasonable" choices of Hi, and for "most 
instances" Hf belonging to hard optimization problem classes, this annealing 
procedure leads to a quantum speedup with respect to classical algorithms. A 
precise meaning can be given to the expression "most instances" by consider- 
ing ensemble of random instances. Continuing with the example of the graph 
coloring problem defined above, one can for instance define probability laws 
on the set of all graphs of N vertices, the most famous one being the Erdos- 
Renyi random graph [30] in which the M edges are chosen uniformly at random. 
Then a property holds for "most instances" if its probability with respect to the 
choice of the random graph goes to 1 in the large size (thermodynamic) limit 
N — > oo. Such ensembles of random optimization problems were actually in- 
troduced in computer science [31] as generators of hard problems on which to 
benchmark classical algorithms. Since then an intense research effort was de- 
voted to their study, in theoretical computer science and discrete mathematics 
of course, but also in statistical mechanics. Random optimization problems can 
indeed be handled by methods first devised for the study of disordered physi- 
cal systems, spin glasses in particular [32]: renaming energy the cost function, 
optimization amounts to low temperature statistical mechanics (one can view 
for instance the graph coloring problem as an antiferromagnctic g-states Potts 
model), an optimal configuration becomes a ground state, and the randomness 
in the instance corresponds to the quenched disorder of spin glasses. This in- 
terdisciplinary approach turned out to be very fruitful, and in the last decade a 
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detailed understanding of the shape of the configuration space of random opti- 
mization problems was reached thanks to the non-rigorous methods of statistical 
mechanics (some of these predictions were later on put on a rigorous mathemat- 
ical basis). In the thermodynamic limit these problems undergo several phase 
transitions when some control parameter of the random ensemble (for instance 
the finite ratio M/N in the case of the graph coloring) is varied, in particular 
the set of ground states gets split into a large number of clusters of close-by 
configurations, the clusters being well separated one from the other in the con- 
figuration space. This understanding also allowed to devise specific ensembles 
of random instances where one can "hide" an arbitrarily chosen unique ground 
state, that remains hard to find if no direct information is available on the hid- 
den configuration. This is particularly useful in the context of the quantum 
adiabatic algorithm, which has most often been studied on instances with a 
Unique Satisfying Assignment (USA). 

As explained above the random optimization problems provided useful bench- 
marks for classical algorithms, it is thus natural to test the efficiency of the quan- 
tum adiabatic algorithm on them, and indeed one of the first proposals [20] stud- 
ied such random instances. Unfortunately simulating quantum computers on 
classical ones is a hard computational task because the dimension of the Hilbert 
space grows exponentially with the number of qubits, hence the numerical inte- 
gration of Schrodinger equation, or the exact diagonalization of the time- varying 
Hamiltonian is restricted to rather small system sizes, whereas computational 
complexity theory classifies the difficulty of problems in the infinite size limit. 
However random optimization problems, viewed from the perspective of sta- 
tistical mechanics of disordered systems, are mean field systems (there is no 
finite-dimensional lattice underlying their definitions) and are as such amenable 
to an analytic resolution. It is thus possible to build upon their classical statis- 
tical mechanics studies in order to include the quantum effects induced by the 
interpolation procedure at the core of the quantum annealing procedure. In par- 
ticular one can investigate the fate of the classical phase transitions mentioned 
above when quantum effects arc added; when these become quantum phase tran- 
sitions [33] as a function of the time parameter of the adiabatic interpolation, 
energy gaps close in the thermodynamic limit and this sets a lower bound on 
the running time of the algorithm, as the adiabatic criterion has to be fulfilled. 
It is thus of a crucial importance to understand the quantum phase transitions 
of random constraint satisfaction problems in presence of quantum fluctuations, 
and in particular their order: generically second order phase transitions are 
associated to polynomially small gaps, while first order transitions (which are 
commonly found in quantum mean field spin glasses [34, 35, 36, 37, 38]) cause 
exponentially small gaps, and in consequence an exponentially long evolution 
time is required for the adiabatic criterion to hold. Another, distinct, mech- 
anism for the appearance of small gaps was pointed out in [39, 40], based on 
an analysis of the perturbative effects of quantum fluctuations on the classical 
energy levels of optimization problems. Some variants of the quantum adia- 
batic algorithm were claimed to circumvent the effects of these "perturbative 
crossings" in [41, 42, 43]. 
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In this brief presentation, as in most of the literature, the quantum adiabatic 
algorithm is viewed as an algorithm to find the ground state of an Hamiltonian, 
or in computer science terms to solve exactly an optimization problem. One can 
however think of it more generically as an approximation algorithm [44]: if the 
allowed evolution time is smaller than required by the adiabaticity criterion then 
the system ends up in an excited state, corresponding to energies higher than 
the global minimum of the Hamiltonian (this is, in physical terms, reminiscent 
of the Kibble-Zurek problem, see [45] for a recent review). But it might be 
that a good compromise can be found between short execution times on the 
one hand, and small excitation energies on the other hand. This would be as 
important from a complexity point of view as being able to find the exact ground 
state. Indeed for several optimization problems it is computationally hard to 
find an approximate value of the ground state energy, and for some of them it is 
hard even to make an estimate better than the energy of a configuration chosen 
uniformly at random [46]. 

In this paper we shall review and extend recent works on the behaviour 
of the quantum adiabatic algorithm on random optimization problems. As we 
explained above this is strongly related to the understanding of the low tem- 
perature phases of quantum spin glasses. The paper is organized as follows. In 
Sec. 2 we shall make a brief introduction to classical and quantum computational 
complexity theory, and define more precisely the quantum adiabatic algorithm. 
Sec. 3 contains a review on classical random optimization problems, their phase 
transitions and their relations to spin glasses. Special attention will be given in 
Sec. 3.9 to the problem of generating random instances with prescribed prop- 
erties, in particular to ensure the non-degeneracy of the ground state (USA 
instance). In Sec. 4 we discuss the thermodynamic properties of quantum spin 
glasses, concentrating in particular on their low energy properties and quantum 
phase transitions, without entering into technical details. The latter are touched 
upon in Sec. 5, where we present several methods, both analytical and numer- 
ical, for the study of quantum disordered systems. Some of these methods are 
then applied to a few representative examples of random optimization problems 
subject to quantum fluctuations in Sec. 6. We finally draw our conclusions in 
Sec. 7. 

In addition to its review character this paper contains original material: in 
Sec. 4.4 and 6.2 we present some details and additional results of two works that 
previously appeared as letters [47, 48]; among the results of these sections the 
study of the gap in presence of an exponential degeneracy of the ground state 
in Sec. 6.2.2 should have a general relevance. The discussion of the quantum q- 
coloring (or antifcrromagnctic Potts model) in Sec. 6.3 was not published before 
and will be further developed in a forthcoming publication [49]. In Sec. 5.6 
we propose a method to extract the gap from Quantum Monte Carlo (QMC) 
numerical simulations, that, to the best of our knowledge, was not discussed 
previously. The discussion on the generation of USA instances of Sec. 3.9 also 
bears some originality in the quantum context. Finally, in Sec. 5.6.3 and 6.2.3 
we show how one can use QMC simulations to detect the clustering transition 
(to be introduced in Sec. 3) of quantum models. 
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Despite its length this review has no pretension of exhaustivity; complemen- 
tary point of views on the quantum adiabatic algorithm can be found in the 
reviews [21, 22, 23, 24, 50, 51, 52] and references therein. 



2. Classical and quantum computations 

2.1. Classical computation theory 

2.1.1. Examples of optimization problems 

We shall give in this section a brief introduction to the classical theory of 
computational complexity [1, 3, 2], and set up some notations that we shall 
use in the rest of the paper. For concreteness we will concentrate on com- 
putational tasks related to combinatorial optimization problems. We shall 
thus consider a discrete configuration space of N variables denoted o~\ , . . . , un , 
each of them taking values in a finite set x, and denote a global configu- 
ration a = (ai, ... ,ctn) € x ■ In most computer science applications the 
variables considered arc boolean, and one usually takes x = {True, False} or 
X = {0, 1}. For consistency with the conventions in vigor in physics we shall 
also use x = { +1 , — 1} , the translations between the various conventions being 
straightforward. The computational tasks we are interested in are defined in 
terms of a cost function that assigns to each configuration g_ a real number. We 
will call this cost the energy of the configuration, or the value of its Hamiltonian, 
and denote it E(a). Let us give some examples: 

• The graph g-coloring problem, in short g-COL, was mentioned in the in- 
troduction and is formalized as follows. Given a graph G = (V, L) with 
V a set of N vertices and L a set of M edges between pairs of vertices, 
one takes x = {lj • ■ • ?<?}: with q > 2 an arbitrary integer, so that each 
configuration a corresponds to the coloring where vertex i is given the 
color Oi. The cost function is 



where the sum runs over all edges of the graph, and S denotes the Kro- 
necker symbol. The cost function thus counts the number of monochro- 
matic edges in the configuration a. 

The following examples involve binary variables that, as explained above, we 
encode with Ising spins, x = {+!> — !}■ 

• The /c-XORSAT problem is defined on a fc-hypergraph G = (V,L): each 
of the M hyper-edges L involves a fc-uplet of variables, with k > 2, thus 
generalizing the notion of usual graphs that corresponds to k = 2. We 
label the hyper-edges with an index a = 1, . . . , M, and denote i\, . . . , %\ 
the indices of the vertices linked by the a-th hyper-edge. In addition to 
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(1) 



the hyper-graph the problem is defined by M constants J a € {+1,-1}, 
and the cost function reads 



fc 

M 1 - J a II a H 

E(<J) = J2 • (2) 

a=l 

It is easily seen that this sum equals the number of hyper-edges a for which 
the condition Cji ... = J a is violated. There are various equivalent 
interpretations of this condition; in the language of coding theory [53] 
this is a parity check rule. By associating Ising spins to {0, 1} variables 
according to at = (— l) Xi it is also equivalent to a linear equation of 
the form x^i + ■ ■ ■ + = y a , where J a = ( — l) Va and the additions 
are interpreted modulo 2. Finally it can be seen as a condition on the 
eXclusive OR of k boolean variables {True, False}, hence the name of the 
problem. 

In the fc-SAT problem one is given an hypergraph and, for each hyper- 
edge, k constants J^,...,J^ g {+1,-1}; the cost of a configuration is 
then defined as 

M k 1 _ tj" 

^sn^Hr*- (3) 

a=lj=l 

Each term of the sum is equal to 1 if, for all the k vertices involved in 
that hyper-edge, one has cr-j ^ J° a ; on the contrary it vanishes as soon as 
one of the k vertices fulfill = J 3 a . In terms of boolean variables this is 
the disjunction (logical OR) of k literals, that are equal to a variable or 
its logical negation depending on the sign of J J a . 

Another example is the so-called l-in-3 SAT (or Exact Cover) problem, 
defined on a hypergraph of triplet of vertices with the cost function 



E(2) = £ 



^ 5 — (Tjl — (7.2 — (T;3 + <7jl(7j2 + Gi\Ui3 + OilOil + 3(7jl <7,-2 (7-3 



a=l 

_(4) 

Each term of the sum is equal to or 1, the former case being realized 
if exactly one out of the three variables involved is equal to -1, the two 
others being equal to 1. 

Note that in all the examples above the energy function is constructed as a 
sum of M indicator functions that take the value (resp. 1) if some constraint 
involving k variables is satisfied (resp. unsatisfied). These examples thus belong 
to the class of Constraint Satisfaction Problems (CSP). In this context one often 
calls a clause each of the individual constraint, and formula the conjunction of 
all the constraints. A formula is said to be satisfied by an assignment a of the 
variables if and only if all the individual constraints are satisfied. A formula is 
satisfiable if and only if there exists at least one configuration that satisfies it; 
in physical terms this correspond to the ground state energy being equal to 0. 



2.1.2. Classical complexity classes 

Given an arbitrary cost function E(a) on a discrete configuration space, one 
can define various computational tasks: 

• The decision task is to answer yes or no to the question "is there a config- 
uration a whose cost is smaller or equal than a given constant C?" . In the 
context of CSP one can further specialize this question by taking C = 0; 
the question thus becomes "is the formula satisfiable?" 

• The optimization task is to compute the minimal value of E over the 
configuration space; the output is thus a real number instead of the yes/no 
answer of the decision task. 

• One can also ask to compute the number of configurations of minimal 
energy (a counting task). 

• Another task is to output explicitly one configuration of minimal energy; 
this could either be any such configurations, or one could require in ad- 
dition that the output configuration is a random configuration, with for 
instance the uniform distribution over all configurations of minimal energy 
(this is a sampling task). 

The goal of computational complexity theory is to classify the difficulty of 
these tasks, in terms of the time and space (memory) requirements of algorithms 
that perform them. Let us emphasize some subtleties in the vocabulary to be 
used: a problem, is, in its loose sense, a set of cost functions. For instance 
the g-coloring problem means all the functions E(a_) defined in Eq. (1), for all 
possible graphs. To be more precise one has to indicate, along with the set of 
cost functions, the version of the problem, among the decision, optimization, 
and counting variants defined above. Finally an instance of a problem means 
one representant of the class of cost functions it includes. An instance of the 
g-coloring problem is thus defined by a graph. 

Let us concentrate first on the decision problems. The NP (standing for Non- 
deterministic Polynomial) complexity class contains the problems for which it 
is easy, for every instance, to check the correctness of the "yes" answer, if the 
algorithm provides as a certificate a configuration er with E(g_) < C. In other 
words for NP problems computing the value of E(a), given a, is by itself an easy 
task, which means a task that can be performed with a number of operations 
growing only polynomially in the size of the input. This is indeed the case for 
all the examples we have given above (the size of the input being here controlled 
by N and M). Of course, even if the answer is easy to check a posteriori, this 
does not mean that the certificate is easy to find a priori. This is true only for 
a subset of the problems in NP, the so-called P (for Polynomial) problems. An 
example of a problem in P is deciding the satisfiability of XORSAT formulas: 
thanks to the mapping onto a problem of linear equations, for any choice of 
the hyper-graph and of the constants J a one can use Gaussian elimination and 
check in a number of operations growing as N 3 whether or not there exists a 



10 



configuration satisfying all the constraints, thus answering the decision question 
with C = 0. The decision versions of 2-SAT and 2-COL, with C = 0, are also 
in P. On the contrary for fc-SAT or g-COL with k, q > 3, no algorithm is known 
to answer the decision question for all possible instances in polynomial time 
(this can of course be done in an exponential time just by inspecting all the 2 W 
configurations one by one) . In fact it is strongly believed that no such algorithm 
can exist: fc-SAT and g-COL (with k, q > 3) belong to the NP-complete subset 
of problems in NP, that are the hardest problems of NP in the sense that any 
instance of any problem in NP can be translated (with a polynomial overhead) 
to an instance of a NP-complete problem. Exhibiting a polynomial algorithm 
for a single NP-complete problem would imply that P=NP, such a collapse of 
the complexity class being held as rather improbable. 

The NP class and its subsets we have briefly discussed is only one example 
among a very large number of complexity classes; let us emphasize that NP 
concerns only the decision tasks whose answer is a yes or no. In general the other 
versions of the problem, which must return a number or a configuration, may fall 
into other complexity classes. In some cases the optimization can be essentially 
reduced to the decision version: if the possible costs are bounded and discrete, 
one can simply make a dichotomy on the value of C and finds the optimal value 
of the cost function by calling the decision problem a number of times growing 
logarithmically with the number of possible values of the cost. But counting 
problems are not reducible in this way, and belong to other complexity classes, 
known as #P and its variants. 

One should also keep in mind that the easiness of a class of problems for 
the decision task does not imply that the other tasks, on the same problem, are 
easy as well. The XORSAT problem is very illuminating in this respect: even 
if its decision version is in P when C = 0, thanks to the Gaussian elimination 
algorithm, solving the decision problem with C > or finding the optimal cost 
of an arbitrary instance takes in the worst-case an exponential time (in the sense 
that no polynomial algorithm is known that performs this task). Indeed, if the 
Gaussian elimination shows that there are no configurations satisfying simulta- 
neously all the constraints, it gives no clue on how to find optimal assignments 
of the variables. The situation is actually even worse: not only it is hard to find 
exactly the optimal cost, but even finding a good approximation for it is also 
hard. This question of approximate resolution of optimization problems is an 
important issue in computer science [44] . In the case of fc-SAT and fc-XORSAT 
it was shown in [46] that finding an approximation of the optimal energy which 
is more accurate than the one obtained by taking uniformly at random a con- 
figuration of the variables is harder than any NP-complete problem. 

Let us finally comment on the notion of execution time of classical algo- 
rithms. In the formal studies of computational complexity theory this time is 
measured in the units of the number of steps performed by a so-called Turing ma- 
chine to execute the algorithm. The Turing machine is a very simplified model 
of a "computing machine" , very far from the complexity of today's computer. 
Its power as a formal tool for the analysis of algorithms relies in its universality, 
formulated in the Church- Turing hypothesis: all classical computers can be em- 
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ulated by a Turing machine with an overhead that grows only polynomially with 
the size of the input to the algorithm, hence the distinction between polynomial 
and exponential run-time is independent of the precise computational model. 

2.2. Quantum computation theory 

We shall now describe a few aspects of the quantum computation theory, and 
contrast it with the classical one. This review being theoretical in nature we 
shall not address the experimental challenges for building quantum computers 
(known as the DiVinccnzo criteria [54]) and only provide a few references to 
experimental works. 

2.2.1. Quantum circuits model and examples of quantum algorithms 

We shall now give a brief presentation of the basics of quantum computer 
science [6, 7, 8], assuming knowledge of the laws and notations of quantum 
mechanics. An introduction to quantum computer science written by and for 
physicists can be found in [55]. The paradigmatic shift from classical computer 
science is the assumption that the elementary variables at the core of the com- 
puter behave quantumly: instead of bits which can take either the value or 
the value 1, one deals with qubits which can be in a coherent superposition of 
the two values. Let us introduce some notations: for a system of N qubits we 
denote T~L the Hilbert space spanned by the ortho normal basis {|er) : g_ G x N }- 
This basis, indexed by the classical configurations, is called the computational 
basis. If x nas d > 2 elements one often speaks of qudits, to emphasize the 
^-dimensionality of the Hilbert space of a single clement. According to the laws 
of quantum mechanics the state of the system is described by a vector \ip) of 
this Hilbert space, i.e. a (complex) linear combination of the vectors of the com- 
putational basis, which has norm 1. The state of the computer evolves during 
the execution of an algorithm; according to the laws of quantum mechanics this 
evolution is represented by the action of a linear operator on the Hilbert space, 
\tjj) — > U\ip), where the linear operator U must be unitary in order to conserve 
the norm of l^). Every quantum algorithm thus corresponds to an unitary op- 
erator; in principle this operator acts on all the qubits of the system, making a 
practical implementation of non-trivial algorithms a seemingly impossible task. 
Fortunately it has been shown [56, 57, 58, 59] that any unitary operator can 
be factorized (with arbitrary precision) as a product of simple operators, called 
gates in this context, that act only on one or two qubits (this is similar, in 
the classical case, to the reducibility of any Boolean function as a combination 
of Not AND gates). Moreover there exist universal sets of gates that contain 
only a finite number of operators. For example, in the case of binary qubits, it 
is enough to take as one-qubit gates the operators P and A, defined by their 
matrix representation 

i.e. P adds a phase of 7r/4 between the two states of the qubit, while the 
Hadamard gate A converts the two vectors of the computational basis into their 
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Figure 1: An example of a quantum circuit on three qubits. In the first step the two first 
qubits are submitted to the phase operator while the third one is acted upon by the Hadamard 
gate. Then <T2 is submitted to a NOT controlled by (73, and finally the system is acted on by 
a (here unspecified) three qubit gate. 



symmetric and antisymmetric linear combinations. The only two-qubit gate that 
completes this universal set is called the controlled NOT (cNOT) gate, that acts 
on the computational basis of the two qubits as \ai, 02) — > ci © 0-2), where 
we use in this section the convention x = {0> 1}> and © denotes the addition 
modulo 2. In other words the cNOT gate leaves the second (controlled) bit 
constant if and only if the first (controlling) bit is 0. 

Quantum algorithms can be conveniently represented graphically by quan- 
tum circuits: the unitary operator U that encodes the algorithm can be factor- 
ized as the consecutive applications of simpler unitary operators, acting possibly 
on a subset of the whole qubits. These elementary operators can be written as 
products of the universal gates displayed above, but this is not compulsory and 
it is often simpler to describe the action of a gate on a large number of qubits 
than its decomposition on one and two qubit gates. An example of quantum 
circuit is shown on Fig. 1. 

Let us now describe some quantum algorithms that exhibit a velocity gain 
with respect to classical computations. The simplest ones shall deal with binary 
functions f(a) from {0, 1}^ to {0, 1} M . In the quantum setting these functions 
are implemented as unitary linear operators [//; note that an unitary transfor- 
mation is invertible, hence Uf must somehow keep trace both of the input g_ and 
of the output /(a) of the function /. A convenient way to fulfill this request is 
to let Uf act on the Hilbert space of N + M qubits, its action being defined on 
the computational basis as 

Uf\g:,a , ) = \a:,a'®f(a)} , (6) 

where © is here the bitwise addition modulo 2. We shall use equivalently the 
notations |ct, <r') and with a = (<ti, . . . , <7at) and a' = (cr^ , . . . , cr' M ), for 

the computational basis vectors, the second notation emphasizing the tensorial 
product between the input and output qubits. At a first (too optimistic) look, 
the laws of quantum mechanics allow to treat in a "parallel" way the 2^ possible 
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inputs of the function /; suppose indeed that the quantum computer is prepared 
in the state 



2 N/2 

where = (0, ... ,0), which can be reached by the application of Hadamard 
gates on the first N qubits to the initial state |0, 0). Then Uf transforms this 
state into 

which seems indeed to contain all the information about the behaviour of / on 
its 2 N possible inputs. However this information is not reachable by an observer, 
because of the measurement axioms of quantum mechanics. A measurement of 
the qubits in the state written above leads to nothing but a random choice 
of a single configuration g_ among the 2 N possible ones, and to the associated 
value f(g_). This trivial observation explains why some thought has to be put 
in devising quantum algorithms that outperform classical ones; the availability 
of linear superpositions is not enough for that, one has to use in a more clever 
way the possibility of interferences between states. 

The simplest example of this strategy is Deutsch's algorithm [5]. Given a 
function / from {0,1} to {0,1} (i.e. N = M = 1), the task is to determine 
whether / is constant or not. Classically one cannot avoid the computation of 
both /(0) and /(l) to answer this question. On a quantum computer this task 
can however be performed with a single application of Uf. Indeed, starting from 
the state |0, 1) and applying Hadamard gates to both qubits leads to 

I(|0) + |1))(|0)-|1». (9) 

Applying the operator Uf, followed by the Hadamard gate on the first qubit, 
produces the state 



|0 J/(Q))- l/(i)) + -l/(Q)) + h J/W)- l/W) + l/(i)) -l/(Q)) 

where we denoted • the logical negation (0 = 1, 1 = 0). If / is constant the 
second term vanishes, otherwise it is the first term that cancels out; measuring 
the state of the first qubit thus yields, without any probability of error, the 
answer to the question. 

In the previous example the "quantum speedup" was rather modest, reduc- 
ing the computational cost from two classical evaluations of / to one application 
of Uf. However, it illustrated the essential ideas behind the much more impres- 
sive gains of quantum algorithms that have been developed later on, and that 
we shall only sketch here. 

One of the problems solved by Deutsch and Jozsa in [9] concerns functions 
/ from {0,l} 7V to{0,l}, with N arbitrary but with the promise that / is either 
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constant or balanced (i.e. takes the value on exactly 2 N 1 distinct inputs). 
Their quantum algorithm (in the refined formulation of [60]) decides between 
these two alternatives with a single application of Uf, and without possibility 
of error, whereas a classical algorithm needs 2 Ar ~ 1 + 1 evaluations of / before 
a definite answer can be given; however a classical randomized algorithm can 
answer after a finite (with respect to N) number of evaluations with a probability 
of error arbitrarily small. 

Simon's algorithm [10] is given as an input a function / from {0, 1}^ to 
{0, 1}^, with the promise that either / is bijective, or that it is "periodic" in 
the (unusual) sense that there exists a binary string t / such that /(a) = 
f{s!) ^ 3- = S.' ©L with again © the bitwise addition modulo 2. This quantum 
algorithm decides between these two alternatives, and allows to determine r in 
the second case, with an expected number of applications of Uf growing linearly 
with N. This has to be contrasted with the exponential number of evaluations 
of / that are necessary for a classical algorithm to solve the same problem. 

A crucial step in Simon's algorithm is a unitary transform known as a Quan- 
tum Fourier Transform. This idea was also exploited by Shor in [11] to devise a 
quantum algorithm for finding the period r of a function / from Z to Z, where 
now the term period is used in a more usual sense: f(x) = f(y) <^> x = y mod r 
(with the promise that r is smaller than some given integer). The importance 
of this result stems from its consequences in the context of number theory, and, 
in a more applied way, to cryptography. As a matter of fact the problem of 
factorizing integers is reducible, via arithmetic theorems, to the period finding 
problem that Shor's algorithm solves in an efficient (polynomial) way. Moreover 
the security of the famous RSA [61] public-key protocol of cryptography is based 
on the inexistence of an efficient classical algorithm for integer factoring. Hence 
the construction of a large quantum computer would have drastic consequences 
for the security of encrypted communications (see [62, 63, 64, 65] for small-scale 
experimental demonstrations) . From the more theoretical point of view of com- 
putational complexity the factoring problem (in its decision version, i.e. given 
N, M two integers, is there p with 1 < p < M such that p divides N, which 
can be used to exhibit a factor of N via a dichotomy on M) is most likely in 
an intermediate difficulty class, namely in NP but outside P and NP-complctc 
(note that the primality decision problem, i.e. the existence of a factor p of N 
was relatively recently shown to be in P [66] , yet without the condition p < M 
this does not solve the factoring problem). Hence the efficient quantum algo- 
rithm devised for solving the factoring problem cannot be used, via reductions, 
to solve all the NP problems. 

Let us also mention another quantum algorithm that is unrelated to those 
mentioned above, and that can be described as follows. Let f T be the function 
from {0,l} Ar to{0,l} that maps all its 2 N inputs to 0, except one fixed string r 
that is mapped to 1. This can be interpreted as an unsorted database with one 
single marked element. In order to discover the value of r a classical algorithm 
(even randomized) cannot do better than computing the value of the function 
on 0(2 N ) inputs. On the contrary Grover exhibited in [12] a quantum algorithm 
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that solves this problem in a number of steps of order 2 N I 2 , i.e. with a quadratic 
speedup with respect to the classical execution time. In fact this problem can 
be reinterpreted as the search for the ground state of the operator Hf = I — 
|r)(r| (where / is the identity operator) on the Hilbert space of N qubits. 
Grover's algorithm works by successive applications of the operators (—1)1^)^1 
and Hi = I — |* )(^'ol (these notations will be useful in Sec. 2.3.5 where we 
shall come back on this problem), where l^o) = 2~ Ar / 2 ^ (T i s the uniform 
superposition of all the states of the Hilbert space. Hi connects any two vectors 
lil)j|<l')j allowing for a "quantum diffusion" between states. The important 
point is that the convergence can be guaranteed within 2 Ar / 2 applications of each 
operator [12], allowing for the quantum quadratic speedup. Grover's algorithm 
is known to be optimal [67] and has been experimentally tested, see [68] and 
references therein. 

Finally, other quantum algorithms have been developed to solve systems of 
linear equations, see in particular [69, 70]. 

2.2.2. Quantum complexity classes 

The classification of problems according to their computational complex- 
ity presented in Sec. 2.1.2 relied on the Church- Turing hypothesis, namely the 
equivalence (within polynomial reductions) of all classical computing devices. 
We shall now briefly sketch the analogous classification that has been devel- 
oped [13, 14, 15], taking as a computing model a quantum computer operating 
algorithms described by quantum circuits. 

The class BQP contains the problems that can be solved with a quantum 
circuit containing a polynomial number of gates; for instance the existence of 
Shor's algorithm [11] demonstrates that the factoring problem belongs to the 
BQP class. This is the quantum analog of the P class, or more precisely of the 
BPP class; indeed the measurement process at the end of a quantum computa- 
tion induces in general some probability of error, that is required to be Bounded 
with respect to the size of the input in the BQP class. 

The closest quantum analog of NP is known as the Quantum Merlin Arthur 
(QMA) class of problems. Let us recall that the (rough) definition of NP we 
gave was the class of decision problems for which a yes answer has certificates 
that can be efficiently checked; for the examples of Sec. 2.1.1 a certificate could 
be provided by a classical configuration a, for which the computation of the en- 
ergy E(q_) was an easy task. In an interactive definition Merlin is the provider 
of the answer and its certificate (using for instance a non-deterministic Turing 
machine) while Arthur is the checker of the certificate. In the quantum transpo- 
sition of this definition Merlin is allowed to give as a certificate of his answer an 
element of the Hilbert space of a quantum computer with a polynomial number 
of qudits, and Arthur can apply a quantum circuit with polynomially many gates 
to this vector in order to verify its validity. Because of the inherent stochasticity 
in the quantum measurement processes some error tolerance has to be included 
in the precise definition of QMA [15], namely the yes instances must have at 
least one certificate that will be accepted by Arthur with a probability close to 
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1, while for the no instances Arthur should be able to reject all the certificates 
that Merlin could try with again a probability close to 1. 

In computational complexity theory the notion of completeness plays a cen- 
tral role: for instance the NP-completc problems are the hardest of the NP 
ones, and as such contains the essence of the difficulty of this class. Similarly in 
the quantum context the QMA-complete problems arc those problems to which 
any member of QMA can be reduced (within a polynomial overhead). In order 
to describe some of the known QMA-complete problems, let us first define the 
notion of local Hamiltonians. Consider the Hilbert space of N qudits spanned 
by {|ct) : a <E x N } (with \x\ finite but possibly > 2). A fc-local Hamiltonian is a 
Hermitian operator H acting on this space, that can be written as H = X) a =i ^ a 
where each H a acts on at most k qudits among the N. The decision problem 
associated to H is to determine whether its smallest eigenvalue (ground state 
energy) is cither < a or > b, where a < b are two given reals and with the 
promise that one of the two alternatives is true (i.e. the smallest eigenvalue of 
H is not in the interval [a, b]). The first QMA-completeness result was obtained 
in [15], where it was shown that the 5-local Hamiltonian problem is indeed 
QMA-complete (for simplicity we keep understood some necessary hypothesis 
on the norms of the H a and on the size of the promise gap b—a). This first result 
was then strengthened in a series of works, that showed the completeness of the 
3-local Hamiltonian problem [71] , then of the 2-local Hamiltonian problem [72] , 
and finally of the 2-local Hamiltonian problem with the further restriction that 
the local interactions H a only couple nearest neighbor qudits on an unidimen- 
sional lattice [73]; this last result only holds if the internal dimension \x\ of the 
qudits is at least 12. 

The locality condition is reminiscent of the form of the classical cost func- 
tions (1), (2), (3), (4) for CSP, that also takes the form of a sum of terms 
acting on a small subset of variables. Consider in particular the fc-SAT prob- 
lem defined in Eq. (3): this is precisely a fc-local Hamiltonian, diagonal in the 
computational basis of N qubits, with each term H a a projector onto the state 
(— J„, • ■ ■ , — J a) of the k qubits This observation triggered the study 

of a quantum generalization of the fc-SAT problem, known as fc-QSAT, where 
the H a are arbitrary projectors acting on k qubits. This problem was first in- 
troduced in [74], where it was shown that the case k = 2 is easy (even on a 
classical computer) while for k > 4 it falls in the QMA-complete class. Several 
results on this problem, and in particular its random version, can be found in 
the original papers [74, 75, 76, 77, 78] and are reviewed in [55]. 

2.2.3. Other approaches to quantum computation 

In the above presentation we have described quantum computations in terms 
of quantum circuits, i.e. the successive action of unitary operators (gates) acting 
on a few qubits. Several alternatives to the quantum circuit strategy have 
been proposed, and the main focus of this review is one of them, the quantum 
adiabatic algorithm. Before presenting it in more details let us just name a few 
of the other perspectives on quantum computation, some of which having been 
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proven to be as universal as the quantum circuit model. 

In the topological quantum computation scheme (see [79] for a review) the 
quantum states that shall be used as qubits are encoded in non-local (topologi- 
cal) degrees of freedom, that increase their tolerance to local decoherence caused 
by an imperfect decoupling from the environment. Experimental realizations of 
this scheme have been proposed, exploiting the non-abelian anyonic statistics 
of excitations in fractional quantum Hall states. 

Another proposal, named "one-way quantum computation" [80], relies cru- 
cially on two exquisitely quantum properties, i.e. entanglement and projective 
measurement. In this scheme the system is initially prepared in an highly en- 
tangled state, and this entanglement is used as a resource for computation, that 
proceeds by a succession of projective measurements on subparts of the system. 

Let us finally mention the quantum walk strategy (see [81, 82, 83, 84, 85] and 
references therein) that promotes the classical random walk procedure to explore 
some configuration space to the quantum level, allowing to exploit interference 
effects between the paths followed by the walk. 

2.3. Quantum annealing, or quantum adiabatic algorithm 
2.3.1. Definitions 

In contrast with the generic quantum computation considerations presented 
above, the main focus of this review is a specific quantum algorithm to solve 
optimization problems, namely the quantum annealing or quantum adiabatic 
algorithm. Let us first emphasize that optimization problems arc intimately 
related to low temperature statistical mechanics. Considering the cost function 
E(a) as an energy, the Gibbs-Boltzmann probability law at inverse temperature 
/3 reads 

M-H^, z(ffl=Ee- w . do) 

where the partition function Z(j3) ensures the normalization of the probability 
law. The latter concentrates on the minima of E(a) in the zero temperature 
limit (/3 — > oo). One can set up a short dictionary translating between the 
optimization and the statistical physics vocabulary: 



Optimization 


Statistical Physics 


cost function 
optimal configuration 
minimal cost 
boolean variables 


energy or Hamiltonian 
ground state 
ground state energy 
spins 



In the classical setting this analogy suggested the so-called simulated anneal- 
ing algorithm [86]: in order to find the minima of the cost function E one can 
perform a random walk in the configuration space, with transition probabilities 
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Figure 2: Schematic picture of the thermal and quantum annealing processes. 



respecting the detailed balance condition (reversibility in the mathematical lan- 
guage) with respect to the Gibbs-Boltzmann distribution, with a time-varying 
temperature that is slowly decreased towards zero. If this decrease is slow 
enough thermal equilibrium is ensured at all times, and at the end of the an- 
nealing the system is found in one of the minima of E: thermal fluctuations 
allow to explore the configuration space and to overcome energy barriers be- 
tween local minima. 

The quantum annealing, or Quantum Adiabatic Algorithm (QAA) [16, 17, 
18, 19, 20], exploits a similar idea but with quantum fluctuations (and barrier 
penetration via tunnel effect) replacing thermal ones (see Fig. 2 for a schematic 
representation of this idea) . To define it more precisely let us introduce the oper- 
ator Hf , which acts in the Hilbcrt space spanned by the classical configurations 
: £L €E X N }- F° r anv cost- function E wc define the associated operator Hf, 
diagonal in the computational basis, with Hf\a) = E(a)\a). The state \tp{t)) of 
the quantum computer evolves according to Schrodinger equation, 



(ii) 



where we used a system of units where Planck's constant h is equal to 1, and 
H{t) is the time- varying Hamiltonian of the system. The algorithm shall be run 
during an interval T of physical time, it will thus be more convenient in the 
following to trade the time t with a reduced time s = t/T G [0, 1]. To perform 
a quantum annealing one has to choose another operator Hi and control the 
system in order to implement an interpolation between the initial and final 
Hamiltonians Hi and Hf, for instance linearly (more general interpolations will 
be discussed below). In this way the state of the system evolves according to 



~\4>(8))=H(8)M8)) , 



H(s) = (l-s)Hi + sHf 



(12) 
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If the initial condition | -0 (0)) is the ground state of H[ = H(0), and if T is 
sufficiently large for the adiabatic condition [25] to hold, then for all s the state 
\ip(s)) is close to the instantaneous ground state of H(s). In particular at the 
end of the annealing \ip(l)) is nearly the ground state of Hf , and a measure of 
the N qubits returns an optimal configuration for the cost function E(a). 

This definition leaves a large variety of possible implementations of the quan- 
tum adiabatic idea. In particular the initial Hamiltonian Hi can be chosen 
arbitrarily a priori, with a few conditions: 

• it should not commute with Hf, otherwise the dynamics is trivial. 

• its ground state should be easy to prepare. 



its construction should not rely on a detailed knowledge of the ground 
state of Hf, that is precisely the problem one tries to solve. 



When the classical variables are Ising spins (\ = {+1,-1}) the computational 
basis can be viewed as the basis of common eigenstates of the Pauli matrices 
along the axis z: o\\a) = o~i\g). In this way Hf is obtained very simply from 
E with the replacement Ui — > of. Then a natural choice for H u that fulfills 
the conditions above, is the action of a transverse field in one direction per- 
pendicular to z, say x for instance: Hj = — Y^iLi&f- Let us recall that of 
acts on the computational basis by flipping the i-th spin, af \a) = with 
ay> = (oi, . . . , <7j_i, — <T.j, . . . , ctjv). Note also that if one is interested in the 
decision problem of the existence of zero energy ground states, then one can also 
consider different cost functions that vanish on the same set of configurations. 
This allows to change Hf in order to improve the efficiency of the algorithm, see 
e.g. [42]. 

An experimental realization of quantum annealing, and a comparison of 
its efficiency with respect to thermal annealing, for a disordered Ising system 
in a transverse field can be found in [19]. This study was performed on a 
macroscopic sample that allowed little control on the final Hamiltonian Hf. 
The experiment of [87] concerned a 3 qubit NMR implementation of a quantum 
adiabatic algorithm; more recently [88] claimed to have controlled an 83 qubit 
quantum computer based on superconducting loops. 

2.3.2. The adiabatic condition 

The first appearance of an adiabatic theorem in the context of quantum me- 
chanics can be traced back to early works of Born and Fock [89] , later rephrased 
in more mathematical terms in [90] (see [91, 92] and references therein for more 
recent discussions); its common formulation in [25] states that, for a system 
evolving according to the time-dependent Schrodinger equation (12), in the ab- 
sence of eigenvalue crossings, the system will follow the instantaneous ground 
state in the limit where the total evolution time T tends to infinity. A more 
precise condition can be found in [25]: let us define A(s) = E\{s) — Eq[s) the 
instantaneous gap of the interpolating Hamiltonian that governs the annealing, 

and b(s) = (l ^1 \ which, once divided by A(s), gives the instantaneous 
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angular speed of the ground state's eigenvector relatively to the first-excited 
state's eigenvector. Then the condition 

t> 7aW VsG[M (13) 

ensures that the probability of not finding the system in the ground state of 
H(l) = Hf at the end of the evolution will be of order at most e 2 . We will refer 
to an evolution time T satisfying (13) as an adiabatic time. In general, b[s) 
can be thought as half the difference in slopes between the ground state's and 
first excited state's energies, and has no singular scaling with the system size 
N; therefore, denoting A m ; n = min se [ ,i] A(s) the minimum value of the gap 
during the annealing, (13) can be replaced by the simpler condition 

r»0(^VA- 2 n ). (14) 

The time of the protocol is governed by the minimum gap and by its scaling 
with N. 

The condition (14) obviously breaks down for any time T if the gap of the 
Hamiltonian vanishes (at finite N) for some value of the interpolation parame- 
ter, which is not expected to happen for < s < 1 for geometrical reasons [93]: 
the Hamiltonian H(s) can be seen as a map from [0, 1] to a real space of di- 
mension 2 2N , in which the subspaces of operators with degenerate eigenvalues 
are hyperplanes of co-dimension 2. Therefore, in the absence of additional sym- 
metries, no strict level crossings are to be expected and A m j n remains strictly 
positive. 

A more subtle situation is encountered if the ground state of the final Hamil- 
tonian Hf is degenerate. In this case, the vanishing of the gap for s getting close 
to 1 is obviously not relevant for the adiabatic evolution of the system. The 
basic idea would be to modify the formulas for A(s) and b(s) to consider only 
transitions between continuations of the classical ground states and first ex- 
cited state(s). However, to the best of our knowledge, no precise formulation of 
the adiabatic theorem exists in this context. For the case in which the ground 
state of H(s) is degenerate with the same degeneracy for all values of the inter- 
polation parameter s, sufficient and necessary conditions for adiabaticity have 
recently been proposed in [94] , extending the work of Wilczek and Zee [95] . Note 
that classically, for a certain class of NP-complete problem such as fc-SAT, NP- 
completcness remains if one conditions on instances with a unique solution [96] . 
Hence, for a worst-case analysis, it is meaningful to study the behaviour of the 
QAA on these instances with a Unique Satisfying Assignment (USA). How- 
ever, if one is interested in an ensemble of random instances (an average-case 
study), then one should be careful that USA instances may not be typical for 
the problem considered, as will be further discussed in Sec. 3.9. 

Finally, let us note that worst-case bounds building on the adiabatic theorem 
for diluted spin systems, as the one relevant for the optimization problems 
considered hereafter, were obtained in [24], allowing to prove that, as for thermal 
annealing, the time for adiabaticity is never larger than an exponential in the 
system size. 
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Figure 3: Schematic representation of the eigenvalues of (15) as a function of time s. Note 
the avoided level crossing in correspondence of the minimum spectral gap. 

2.3.3. The finite-time Landau-Zener example 

As explained above, in absence of special symmetries in H(s) and when 
the size of the system is finite (even if very large), true level crossings are not 
expected; but levels may still get extremely close, defining the appearance of 
avoided level crossings. It is then very useful to consider the following "reduced" 
Hamiltonian that describes such an avoided crossing (see Fig. 3): 

«»<•>- (V V-»-))- (15) 



The instantaneous gap is A(s) = 1\/ + b 2 (s — s*) 2 , and even when the two 
diagonal ( "unperturbed" ) elements are equal in s = s* the states do not cross 
but are split by a gap A m i n = 2j. The advantage of this simplified formulation 
is that it is exactly solvable in the limit of an evolution going from s = — oo to 
s = oo, the Landau-Zener formula [97, 98] giving the probability P of a diabatic 
transition to an excited state as P = e~ 2T7rj ~/ b , which has for consequence 
the necessary condition for an adiabatic process T 3> &7~ 2 — 6A~ 2 n , which 
is precisely (13). For evolutions of finite duration, as the ones relevant in our 
context, this formula has to be corrected [99, 100] but the conclusions remain 
unchanged. 

Finally, note that it is possible to extend this formula to consider several 
level crossings [101], and to build on these exact results for this simplified model 
to make predictions for realistic systems involving an extensive or exponential 
number of levels [102, 103]. 

2.3.4- Universality of the quantum adiabatic algorithm 

It could seem at first sight that the quantum adiabatic algorithm has little 
to do with the algorithms based on the quantum circuits model described in 
Sec. 2.2.1; the former is based on a continuous time evolution of the quantum 
computer, and is aimed at finding the ground state of the Hamiltonian Hf, 
while the latter class of algorithms proceed via a discrete succession of unitary 
transformations, and encompass a large variety of computational tasks. An 
equivalence between the two paradigms has however been demonstrated, in the 
following sense. On the one hand, a continuous time annealing procedure can 
be approximated, with an arbitrary precision and with a polynomial overhead, 
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by a series of discrete transformations [93, 26]. In the reverse direction, it was 
shown in [104] that any quantum circuit model can be converted into a quantum 
annealing procedure, by using a final Hamiltonian Hi introduced in [15] whose 
ground state has a positive overlap with the final state of the original quantum 
circuit. Moreover the minimal gap along this interpolation was proven in [104] 
to be only polynomially small in the number of gates of the circuit, hence 
the requested time for the adiabatic algorithm is polynomial in the size of the 
circuit. This result was strengthened in [73], which demonstrated that the 
annealing Hamiltonian can be written with nearest-neighbor interactions on an 
unidimcnsional lattice (the price to be paid being the internal dimension of the 
qudits that has to be larger than 9). 

2. 3. 5. Deficiencies and improvements of the quantum adiabatic algorithm 

The adiabatic theorem stated above provides a very simple and generic con- 
dition under which the quantum adiabatic algorithm is guaranteed to find a 
solution, if any, or at least a minimal energy configuration, to any given op- 
timization problem. However, this does not mean that the algorithm will be 
efficient in finding this answer; in fact its performance will strongly depend on 
the possibly fast closing of the gap along the annealing path. A detailed discus- 
sion of these phenomena is the main focus of this paper. However, it is useful 
to give here a short account of the main points to be discussed in the following. 
As a first example, specific instances of fc-SAT on which the QAA is inefficient 
because of an exponentially small minimum gap were constructed in [26, 27]. 
More generally, we shall see in the following that for random optimization prob- 
lems, two quite general mechanisms may cause gaps closing exponentially fast 
with the system size, and hamper the performances of the quantum adiabatic 
algorithm: 

• The low energy states of the adiabatic Hamiltonian for s = and s = 1 
are very far away one from each other in the Hilbert space. One may 
in particular expect a spin glass phase for s close to 1, and a quantum 
paramagnetic phase for s close to 0, separated by a quantum phase tran- 
sition. Such a phase transition gcnerically leads to a vanishing gap in the 
thermodynamic limit, with a scaling in the system size that depends on 
the order of the transition: in general, the gap closes polynomially fast 
if the transition is second-order, and exponentially fast if it is first order. 
We will come back on such quantum phase transitions in Sec. 4, and in 
particular on their effects on typical constraint satisfaction problems in 
Sec. 6. 

• According to classical spin glass theory (Sec. 3), typical difficult problem 
cost functions are characterized by the existence of many very "different" 
minima (local or global), leading to a very complicated structure for the 
low energy phase of the final Hamiltonian. It may happen that the addi- 
tion of quantum fluctuations leads to many exponentially small gaps be- 
tween different states even within the spin glass phase. This phenomenon 
will be further discussed in Sec. 4.3.2, and on a particular example in 6.3. 
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Before entering the more detailed discussion of the problematics related to 
the efficiency of QAA for solving real optimization problems, let us note that the 
QAA setting is fairly general and leaves open a lot of directions for improvement. 
For instance, we assumed here that H (s) interpolates linearly in time between 
Hi and Hf] but more general interpolation rates can be considered, as will be 
discussed below. Another important freedom in the setting of the QAA is the 
choice of the initial Hamiltonian Hi (and of the hnal one, Hi, if one is only 
interested in the satisfiability decision problem) . The most general formulation 
of a QAA should thus be that of a smooth mapping from some interpolation 
range, that can be taken without restriction to be [0, 1], to the space of Hcrmitian 
operators on (some) Hilbert space, with the constraint that s = 1 is mapped on 
the classical energy cost function, and s = to some operator defined without 
using information on the ground state of Hf, and such that its ground state is 
easy to prepare. Although we do not intend to give a more precise definition 
of these conditions here, their meaning should be clear on concrete examples. 
The latter formulation naturally maps the question of finding the best annealing 
path to a geometrical problem in a Hilbert space [105]. 

Let us now come back to the unsorted database search introduced in 2.2.1, 
to show on this simple example how a modification in the interpolation rate 
can lead to important changes of the adiabatic time of the QAA with fixed 
initial and final Hamiltonians. We will follow here the works of [106, 107, 26]. 
We recall that this problem can be seen as the search for the ground state of 
the classical Hamiltonian H{ = I — |r)(r|, where |t) is some fixed vector of 
the Hilbert space. Let us take for initial Hamiltonian Hi = I — \'&o)( 1 $>o\, with 
l^o) = 2~ Ar / 2 X) CT 1^) the uniform superposition of all the states of the Hilbert 
space defined in 2.2.1. It can be seen that any state \a) — |cr') with a, cr' 7^ r is an 
eigenvector of H (s) with eigenvalue 1 for all s. One can construct 2^ — 2 linearly 
independent such states, that are all orthogonal to both \t) and \*f?o): thus, only 
the subspace spanned by |r) and \^>o) is relevant for the adiabatic evolution. The 
Hamiltonian restricted to this subspace can easily be diagonalizcd, leading to a 
gap A(s) = yfl — 4 (1 — 2~ N ) s(l — s) which is minimal for s = 1/2, resulting 
in A m j n = 2~ N / 2 and in a growth of the adiabatic time proportional to 2^, 
which is also the duration of a naive exhaustive search. However, we know 
from the Grover circuit algorithm [12] that a quantum computer is able to get 
a quadratic speed-up and to find the ground state of Hf in a time growing only 
as 2 JV/ ' 2 . The reason why the QAA seems, in its naive setting, inefficient is 
that the condition (14) is realized only at one particular point of the spectrum 
(s = 1/2) but leads to a constraint on the speed of evolution for all values of 
the interpolation parameter s, even when the gap is large and the annealing 
could be faster without inducing diabatic transitions. Therefore, it is better 
to make a more precise use of the condition (13) and to do the evolution with 
the change of parametrization H(s) = H((p(s)), allowing to vary the speed of 
evolution as a function of the parameter s. Then, with the notations of (13), 
using b(s) = |(l|<iffy<is|0)| = b(ip(s))dip(s) / ds and the simple bound b(s) < 1, 
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the adiabatic condition (13) translates into: 

< eA(^( S )) 2 (16) 

Solving this differential equation as a function of e and N with the boundary 
conditions <p(0) = and y(l) = 1 allows to find the optimal annealing schedule 
<p{s) and fixes the annealing time as T — — , which is the expected quantum 
speed-up. In general it is more difficult to find the best annealing rate for a 
fixed evolution time T\ such questions are related to quantum optimal control, 
as presented in [107, 108, 109]. The intuitive idea is that, at least if the location 
of the gap is exactly known, it should always be possible to trade the scaling 
of the adiabatic time with A~? n of (14) into a scaling with A~[ n , in the same 
fashion as was done for the Grover problem. On the other hand, it is easy to see 
that one cannot do better; in fact, the regime of T <C A~ ■ corresponds to the 
fast passage regime of [25] in which the system is strongly diabatic. In particular, 
one cannot hope to change the scaling of the adiabatic time with the system 
size only by playing on the time-dependence of the evolution Hamiltonian. 

The choice of Hi is more crucial, and for general optimization problems, it is 
mainly an open problem to understand whether the modification of the anneal- 
ing path can change the scaling of the time needed for the adiabatic condition 
to hold. Such a possible change was argued for in [110, 111] to avoid gaps at 
a first order phase transition for fully connected models by introducing a two- 
parameter annealing path (see also [112] for another example of a two-parameter 
annealing path). Alternatively, a randomization of Hi was proposed in [41] to 
avoid gaps of the second type in the classification above, that appear within the 
classical spin glass phase, for a particular problem; but its efficiency for more 
general optimization problems is still an open question. This proposal will be 
further discussed in Sec. 6.2.5. Let us finally emphasize that the existence, for 
a given problem, of an annealing path allowing for a fast adiabatic evolution is 
not enough if the time needed to find this particular path grows exponentially 
fast with the system size [43]. 

2.3.6. Quantum annealing without adiabaticity and approximation issues 

Finally, an important observation is that most of the works up to now focused 
on the efficiency of the quantum adiabatic algorithm in solving exactly the 
problems, that is finding the ground state of the final Hamiltonian Hf . However, 
the question of finding approximate solutions to an optimization problem is of 
great importance, both theoretically [44] and practically. A convenient way to 
quantify the performance of a given algorithm in finding an approximate solution 
to an optimization problem is to introduce its residual energy on a given time 
T, which is the difference between the lowest value of the cost function it can 
achieve in time T and the absolute minimum of the cost function. A zero residual 
energy means that the algorithm can find a solution in time T, while the trivial 
residual energy corresponds to the energy of a randomly chosen configuration, 
that can be achieved with T — 0. Between these two extreme cases, finite (N 
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independent) and extensive (proportional to N) residual energies shall also be 
distinguished. 

This leads to the natural question of how the evolution time must grow with 
the system size for the residual energy to be under a given threshold. Classi- 
cally, it is known that for certain hard problems such as fc-SAT or fc-XORSAT, 
obtaining a non-trivial residual energy can already require an exponentially long 
time [46]. Hence a fast non-adiabatic evolution has a computational interest if 
one can find a good compromise between the evolution time and the residual en- 
ergy. In the classical case, hardness of approximation results can more generally 
be obtained via the PCP theorem [113]. In the quantum complexity literature a 
quantum analog of the PCP theorem has been conjectured in [114]. For recent 
works on the approximation algorithms in the quantum complexity setting we 
refer the readers to [115, 116, 117]. 

Still, the performances of QAA in finding approximate solutions remain 
widely unexplored. Already in [102], it was shown some evidence that QAA 
could outperform classical simulated annealing within the same exponential 
scaling of the running time. To make more general theoretical predictions on 
the residual energy obtained by the QAA, it will be necessary to extend the 
relationship between the spectrum and the behavior of the quantum time evo- 
lution beyond the adiabatic criterion that focuses on the gap above the ground 
state. It has for instance been shown in [103] that the metastable continuation 
of the ground state that emerges after a first order phase transition for fully 
connected mean field models is particularly relevant for quantum evolution on 
sub-exponential time scales, and leads to extensive residual energies for evolu- 
tion times that do not grow as fast as the time for adiabaticity (14). Such a 
connection is expected to have a wider range of validity; in particular to hold 
for random optimization problems as the ones studied hereafter, although no 
quantitative prediction has been obtained yet for these models. 

3. Classical random optimization problems and their connection with 
mean field spin glasses 

3.1. Optimization in the typical case, and the statistical physics of disordered 
systems 

The theory of classical computational complexity [1, 3, 2] that we described 
in Sec. 2.1 considers the difficulty of a problem in the worst-case. For instance 
the fact that g-coloring belongs for q > 3 to the NP-complete class means that at 
present there is no polynomial-time algorithm able to decide the colorability of 
every possible graph. However this does not mean that all the graphs are equally 
difficult, and in fact for many NP-complete problems there exist algorithms that 
do work efficiently on a large set of instances. This raises the question "where 
are the really hard instances of NP problems [118]?" , and how to construct such 
hard instances efficiently 

The idea of using random instances of Constraint Satisfaction Problems 
(CSP) as benchmarks for algorithms emerged in the 80's; however the first 
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ensembles proposed turned out to contain mostly easy instances for known al- 
gorithms, and it was only at the beginning of the 90's that two seminal papers by 
Cheeseman, Kanefsky and Taylor [118] and Mitchell, Selman and Levesque [31] 
introduced the random ensembles that answered positively the question above. 
Instances from these ensembles are actually very simple to describe: in the 
coloring case one creates a graph by selecting uniformly at random M edges 
among the (^) possible ones, i.e. one constructs an Erdos-Rcnyi G(N, M) ran- 
dom graph [30]. The large-size limit (N — > oo) has to be performed with M 
growing like N, in other words the thermodynamic limit for these instances is 
parametrized by a finite real number a = M/N. For fc-SAT random instances 
the construction is generalized to random hyper-graphs, the M fc-uplet of in- 
dices in Eq. (3) being chosen uniformly at random among the (^) possible ones, 
and the signs J 3 a of the corresponding literals are chosen to be ±1 with proba- 
bility 1/2. Again the large-size limit is taken with a = M/N fixed. The authors 
of [118, 31] performed extensive numerical experiments on such randomly gen- 
erated instances. Using complete algorithms they determined the probability 
P(a, N) that an instance with N variables and M = aN constraints has a 
ground state of vanishing energy (i.e. is g-colorable, or satisfiablc depending 
on the case). This probability is obviously a decreasing function of a: it can 
only become harder to satisfy all the constraints as their number is increased. 
What came as a surprise at that time is the fact that for larger values of N 
the probability of satisfiability decreased in a steeper and steeper way, which 
suggested the following satisfiability conjecture: 

1 if Q < Q, , . 

lim P(a,N) = \ , (17) 

N^ca 10 if a > a s 

where a s is some fixed threshold value, that depends on the problem considered 
(coloring or satisfiability), and on the parameters k, q. In more physical terms 
this threshold phenomenon corresponds to a phase transition between a SAT 
(or COL) phase where almost all instances are satisfiable (colorable) and their 
ground state energy is zero to an UNSAT (UNCOL) phase in which almost 
none of them is, and the average ground state energy is positive. Moreover the 
hardest instances, in terms of the time required for the algorithms to decide 
their satisfiability, are those with a a s : for a -C ct s the problem is under- 
constrained, and it is easy to find configurations satisfying all the constraints 
simultaneously, while for a 3> a s there are so many constraints that it becomes 
(relatively) easy again to discover an unavoidable contradiction between them. 

Since their introduction these ensembles have been the subject of a very 
important research effort in computer science, discrete mathematics, and statis- 
tical physics; they have played (and still do) a prominent role in understanding 
the origin of algorithmic hardness. The rigorous works on this problem were 
first aimed at the proof of the satisfiability conjecture (17) and the determina- 
tion of the threshold a s . The main outcomes of this line of research have been 
a proof of a weaker version of (17) where a s is allowed to depend on N [119], 
and rigorous upper and lower bounds on a s [120, 121, 122]. These bounds are 
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asymptotically tight when k,q get large [123]. In addition statistical mechanics 
techniques, starting from [124], have also been applied to these problems and 
have led to quantitative computations of the value of a s [125, 126, 127, 128]. 
Moreover these studies unveiled several new qualitative features besides the sat- 
isfiability transition at the threshold a s ; it has been shown in particular that in 
the SAT phase a < a s there exist further phase transitions [126, 129, 130] that 
affect the organization of the solutions of the random CSP in the configuration 
space, and that are at least as relevant as the SAT-UNSAT transition to un- 
derstand the algorithmic hardness. Note that even if the statistical mechanics 
techniques are not rigorous from a mathematical point of view, many of the 
insights they offered on the features of random CSP have later been turned into 
mathematically rigorous statements [131, 132, 133, 134, 135]. 

Our goal in the remaining of this section is to review the picture of random 
CSP that has been obtained by physics methods (see [136, 137, 138] for textbook 
presentations). Let us first explain in generic terms why statistical mechanics 
is a natural tool for their study, besides the superficial analogy between the 
satisfiability threshold phenomenon and phase transitions of real materials. As 
explained with the dictionary introduced in Sec. 2.3.1, the cost function E(a_) for 
one instance of a CSP can be viewed as an energy function; turning to random 
CSP, this energy function becomes itself a random object. Physical systems 
defined via random constructions have been studied for decades in physics (an 
early example being the Anderson model [139] of localization); in that context 
the randomness in the energy function of one instance (for instance the choice 
of the graph in random coloring) is usually called quenched disorder of that 
sample. Random CSP can thus be studied from the perspective of the statistical 
mechanics of disordered systems. Moreover they belong to the so-called mean 
field class of models, because their structure is unrelated to a finite-dimensional 
physical space: in the Erdos-Rcnyi definition of a random graph all pairs of 
vertices have the same probability to become neighbors (i.e. be linked by an 
edge), there is no a priori Euclidean distance between them. 

In order to make the results on random CSP accessible to readers not ac- 
quainted with the field of statistical mechanics of disordered systems we shall 
make a detour and first discuss simpler models, introducing the necessary in- 
gredients progressively. In Sec. 3.2 we shall introduce the disordered physical 
systems that are most relevant to this discussion, namely spin glasses, and dis- 
cuss the various kinds of mean field models. Then in Sec. 3.3 we present the 
random energy model [140], the simplest disordered model that yet displays a 
phase transition important for the following discussions. In Sec. 3.4 we move 
on to a slightly more complicated model, the so-called fully connected p-spin 
model, and discuss its interpretation in terms of the physics of glasses. We then 
come back to the main focus of our interest, i.e. random CSP; in Sec. 3.5 we 
discuss random instances of the XORSAT problem, followed in Sec. 3.6 by a 
presentation of a toy model that exhibits, in a controlled way, the transitions 
of the random ^-satisfiability and q-coloring model. The latter are discussed 
in Sec. 3.7, without entering into technical details of their derivations, some of 
which will be given in Sec. 5.1. In Sec. 3.8 we will discuss the consequences 
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of these transitions for thermal annealing. Finally in Sec. 3.9 we discuss the 
generation process of random CSP, with a particular interest on ensembles of 
instances with a Unique Satisfying Assignment (USA); these have a special in- 
terest as benchmarks for the quantum adiabatic algorithm. 

3.2. Mean field spin glasses 

Spin glasses can be prepared as alloys of two elements, with a small fraction 
of a magnetic clement (Fe for instance) being added to a metallic host with 
no magnetic properties (Au). This mixture is prepared as a liquid phase at 
high temperature; when the sample is cooled down and becomes a solid the 
position of the magnetic impurities becomes frozen (quenched) to a random 
location. The magnetic moments (spins) carried by the impurities interact with 
one another, but, depending on the distances between them, their pairwise 
interactions can be either ferromagnetic, favoring the alignment of the two spins, 
or antiferromagnetic, forcing them to point in opposite directions. When the 
temperature is varied there appears in these compounds a phase transition for 
the magnetic degrees of freedom. The low temperature phase is an unusual 
state, with frozen moments but no periodic order; hence, the name spin glass, 
by analogy with amorphous window glass, slow to respond to changes in external 
controls, accompanied by non-ergodicity, behaving differently depending on the 
order in which external perturbations, such as magnetic field or temperature, 
are applied. Nowadays, the expression "spin glass" is however used much more 
broadly to refer to systems that exhibit glassiness owing to the combination of 
quenched disorder and frustration. 

In 1975 Edwards and Anderson (EA) introduced in [141] a model for spin 
glass materials, in which the magnetic moments are modeled by Ising variables 
<7j = ±1, lying on a regular finite-dimensional lattice and interacting via the 
energy function 



where the sum runs over the pairs of neighboring spins i and j. The couplings 
Jij are chosen at random from a given distribution (for instance a Gaussian 
one) that allows both positive (ferromagnetic) and negative (antiferromagnetic) 
values for . Statistical models defined on finite-dimensional lattices are very 
difficult to treat analytically, even in the pure case without disorder. The situa- 
tion only becomes worse with the inclusion of disorder, and no analytic solution 
of the Edwards- Anderson model can be hoped for in dimensions larger than 1 . 

The usual prescription of field theory is to start working out the mean field 
version of a model (usually qualitatively correct in large dimensions). For spin 
glasses this was first investigated, after [141], by Sherrington and Kirkpatrick 
(SK) [142]. In the SK model N Ising spins interact with the energy function 




(18) 




(19) 



i<j 
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where the sum is now over all couples i ^ j (making the graph of interaction 
a fully connected, or complete, one). For the thermodynamic limit to be well- 
defined the random couplings Jij must be individually weak, for instance they 
can be chosen to be Gaussian with zero mean and variance of order 1/N. 

This model has played a fundamental role in the theory of spin glasses. 
Despite its mean field character the quenched disorder in its definition makes 
the computation of its free energy a very difficult problem, that was only solved 
in 1980 by Parisi [143], via the development of the replica method in its Replica 
Symmetry Breaking (RSB) form. The SK model exhibits a phase transition 
from a high temperature, paramagnetic phase, to a spin glass phase at low 
temperature, characterized by a proliferation of metastable states in a very 
complex free energy landscape. The reader is referred to the books [32, 144] 
for details on the replica method and the original works on the characterization 
of the spin glass phase of the SK model. The methods originally employed by 
physicists for the resolution of this model were highly non-rigorous. However 
the value of the free energy computed by Parisi was rigorously proven to be 
exact, much more recently, by Talagrand [145], building on the interpolation 
method of Gucrra and Toninclli [146]. 

Let us introduce here some variants of the SK model and set up some ter- 
minology that will often appear in the rest of the discussion. A first twist on 
Eq. (19) consists in promoting the pair-wise interactions to p-wise couplings, 
leading to 

E{o) = - . (2°) 

ii<...<j p 

where the sum is over all p-uplets of spins, and the couplings Jj 1 ...j arc Gaussian 
random variables of zero mean and variance of order l/iV p_1 . This model is 
known as the fully connected p-spin model, and was first introduced and studied 
in [140, 147]; the replica theory developed for the SK model is also applicable 
to this model, that we shall discuss slightly further in Sec. 3.4. The models 
defined in Eqs. (19), (20) have a mean field nature, because each variable at 
interacts (weakly) with all the other variables, destroying completely any notion 
of finite-dimensional distance between the variables; this class of models defined 
on complete graphs are usually called fully connected mean field models. 

But as we already mentioned there exists another class of mean field models, 
dubbed sparse, or diluted, or finitely-connected, in which each degree of freedom 
interacts strongly (i.e. with a coupling of order 1) with a finite number of 
neighbors, the latter being chosen in some random way unrelated to a finite- 
dimensional space. For instance Viana and Bray [148] considered a model of 
pairwise spin glass interactions along the edges of an Erdos-Rcnyi random graphs 
(as defined in Sec. 3.1). More generically a finitely-connected mean field model 
can be defined with other types of sparse random graphs (or hypergraphs to 
include interactions between more than two spins), as long as the connectivity 
(degree) of each vertex remains finite. One way to define a random graph 
probability law is to impose its degree distribution, i.e. to generate uniformly 
at random a graph among those that have a prescribed fraction q$ of isolated 
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vertices, qi of vertices adjacent to a single edge, and so on and so forth. A 
particular case of this construction that will be met often in the following is the 
regular one: a random c-regular graph is a graph chosen uniformly at random 
among all the graphs in which each vertex has exactly c neighbors (or belong 
to exactly c fc-uplets for the hypcrgraph generalization). 

All these sparse random graphs models share a crucial property: they are 
locally tree-like. In other words if one selects an arbitrary vertex i in a sparse 
random graph of size N, with a probability which goes to one in the limit 
N — > oo the shortest loop around i will be larger than any fixed length [30]. 
Statistical mechanics models defined on trees are trivial: they can easily be 
solved by recurrence (somehow like in unidimensional models with the transfer 
matrix method). The richness of the models defined on random graphs comes 
from a subtle combination between their locally tree-like character, and the 
existence of long loops. The latter are very important in creating self-consistent 
boundary conditions, and in avoiding the pathologic surface to volume ratio 
of tree models. Sparse random graphs are sometimes called Bethe lattices, in 
honour of the Bethe approximation that becomes exact on trees; note however 
that this terminology can be misleading, some authors using it as a synonym for 
infinite Cayley trees, some restricting it to the case of regular random graphs. 

From the introduction to random CSP of Sec. 3.1 it should be clear that the 
diluted mean field models will ultimately be more useful in this respect than the 
fully connected ones (though other optimization problems, not described here, 
are defined on complete graphs [149, 150]). They are unfortunately much more 
difficult from a technical point of view. The replica method could be adapted 
to deal with sparse random graphs (see [151] and references therein) but yields 
functional equations under a form that is not directly amenable to numerical 
resolution. An alternative formulation was developed under the name of cavity 
method [152, 153] and allowed to bypass this difficulty. This method, that we 
shall review in Sec. 5.1, yields formally exact predictions for the thermodynamic 
limit of the free energy of models defined on random (hyper)graphs, even if in 
some cases the extraction of actual numbers out of the method can be difficult. 

Before getting to the discussion of the picture of random CSP provided by 
statistical mechanics studies let us discuss, as announced above, simpler models 
of disordered systems. 

3.3. The random energy model 

By definition the energy function E(g_) of a disordered system, by contrast 
with a pure or ordered one, is a random object. For generic local cost func- 
tions (in the sense of Sec. 2.2.2, i.e. that are a sum of terms each involving a 
finite number of spins), the energies of the 2 N configurations (for Ising spins) 
are random variables, correlated one with the other: for instance in the SK 
model (19) the number of independent couplings is only of order N 2 . It is 
however very instructive to study the simplified (but non-local) Random En- 
ergy Model (REM) of Derrida [140], which keeps the random character of the 
energy function but discards the correlations between the energies E(a) of var- 
ious configurations. More precisely, in the REM one assigns to each of the 2 N 
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configurations an energy E(g_) drawn independently at random with a Gaussian 
distribution of density 

P{E) = ^=L . (21) 
V Nit 

The simplest way to solve this model is to use the micro-canonical ensemble. 
Let us denote n(E)dE the number of energy levels belonging to the interval 
(E, E + dE); its average over the realizations of the disorder (the choice of the 
energies) is easily computed: 

n{E) = 2 N P{E) ~ e N{i^E-,N-)) = e Ns(E/N) j (22) 

where ~ denotes here equality at the leading exponential order when N — > oo, 
and the micro-canonical entropy s(e) for the reduced intensive energy e = E/N 
is 

s(e) = log2 - e 2 . (23) 



This function is positive on the interval [eo, — eo], with eo = — -y/log2; for energies 
E corresponding to this interval n[E) is exponentially large and the random 
variable n{E) is thus typically close to its average (with fluctuations of order 

1/2 

n(E) ). On the other hand if E is outside the interval the average number 
n(E) is exponentially small, hence in the vast majority of samples the number 
n(E) is equal to zero. The typical value of the free energy can then be computed 
by the Legendre transform of the typical micro-canonical entropy: 

/RB M = -i lim log / £ ° e N[-Pe+s(e)] de ^ j nf [ e _ Ts(e) ] ? ( 24 ) 
P N ^°° Jeo ee[e ,-e ] 

where we evaluated the integral by the Laplace method. A transition between 
two regimes thus arises at a critical temperature T c such that -r = ds j^ | = 
2-y/log2 and the thermodynamic behavior of the model follows: 

• i) For T < T c , /rem = — Vl°g2 and the system is frozen in its lowest 
energy states (the integral in (24) is dominated by the lower edge eo of 
the integration domain). One can show that only a finite number of con- 
figurations (and only the ground state at T = 0) contribute significantly 
to the partition sum (see for instance [137, 154]). The energy gap between 
them is finite. 



(ii) For T > T c , /rem = — — Tlog2; exponentially many configurations 
contribute to the partition sum. 



The free energy of the model is thus non- analytic at T c . This phase transition is 
often called a "condensation" transition (or Kauzmann transition in the context 
of glasses, see below), because it separates a high temperature phase in which 
the Gibbs-Boltzmann distribution is spread over an exponential number of con- 
figurations from a low temperature phase where this support condenses on a 
much smaller number of configurations. This kind of transition appears, with 
additional subtleties, in many of the more complicated mean field disordered 
systems that we shall discuss in the following. 
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3-4- The fully connected p-spin model 

The next model we would like to discuss is the fully connected p-spin model 
of Eq. (20), studied in particular in [147] (see [155] for an extensive pedagogical 
discussion). We assume p > 3 here, the case p = 2 of the SK model being 
qualitatively different. At variance with the REM, this is a p-local cost function 
(it is a sum of p-spin interactions) and the energies of the various configurations 
are correlated; flipping one of the spins does not completely change the energy, 
some continuity is preserved in the energy landscape of the configurations. There 
is however one common feature with the REM (which is actually the p — > oo 
limit of the p-spin model [140]): it also exhibits a condensation transition at 
some temperature T c , accompanied by a non-analyticity of the free energy. The 
difference is that in the low temperature phase the Gibbs-Boltzmann measure is 
supported by a small number of "pure states" . The latter, that take the place of 
the low energy single configurations of the REM, are whole sets of exponentially 
many correlated configurations. Each pure state therefore has an extensive 
"internal" free energy with both energetic and cntropic contributions. It is not 
easy to give a clear-cut definition of pure states in mean field disordered systems; 
the reader might want to think about them as a generalization of the pure 
states of low temperature ferromagnets with positive/negative magnetizations. 
The partition of the configuration space into pure states has both static and 
dynamic characterizations: long-distance connected correlation functions vanish 
inside one pure state, and the dynamics remains trapped for a long time in the 
pure state it started in. For a given realization of the disorder let us index with 
7 the pure states on which the system is decomposed. The partition function 
can be written as a sum over the pure states, 

z = J2 z 7 , z^y. e ~ PE{2) > a = ~m lo § z 7 . (25) 

7 21^7 

where we denoted / 7 the internal free energy density of the pure state 7. In 
mean field disordered systems, at sufficiently low temperature, there exist ex- 
ponentially many pure states, whose internal free energy density can vary in an 
interval [/ m i n , /max]; one defines a complexity, or configurational entropy, £(/), 
such that e NTi ^ gives, at the leading exponential order, the number of pure 
states with internal free energy /. Then the computation of the free energy is 
a generalization of (24) , 

/ = -~ lim log / /mOX e^-^+^^d/to - inf [/ mt - TS(/ int )] . 

P N^OO Jf m . n /i»te[/mi„,/max] 

(26) 

The condensation transition is thus due here to a competition between the 
internal free energy of the pure states and their degeneracy (configurational 
entropy) ; at low temperatures the integral becomes dominated by the lower edge 
/min of the integration domain, where the configurational entropy generically 
vanishes (note that in general £(/) also depends on the external parameters 
like the temperature), and remains zero at lower temperatures. 
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Another difference between the REM and the p-spin modci is the existence of 
another transition at a higher temperature Td > T c . This so-called "dynamic" 
temperature marks the appearance of pure states inside the Gibbs-Boltzmann 
measure; for higher temperature the space of configurations is essentially con- 
nected and ergodic, only for T < Td the decomposition in pure states is relevant 
and the complexity £(/) non-trivial. This transition has thus a direct impact 
on the dynamics of the system: for T < Td it is not possible for a physical 
dynamics to equilibrate and the ergodicity is broken. However the free energy 
has no singularity at this temperature. 

This model has played a very important role [156, 157, 158, 159] in the de- 
velopment of a first principle theory of the structural glass transition, known 
as the Random First Order Transition (RFOT) theory, see [160, 161, 162, 163] 
for recent reviews. In this scenario, the dynamic transition at Td, with no im- 
pact on the statics, corresponds to the transition of the Mode Coupling Theory 
(MCT) [164], while the condensation transition at T c is an idealization of the 
thermodynamic glass transition envisioned by Kauzmann [165]. 

It is worth mentioning that in many models, a third phenomenon is observed 
as the temperature is further lowered, called the Gardner transition [166]. It is 
a transition towards a more complicated phase, similar to the one found in the 
Sherrington-Kirkpatrick model [32] . 

3.5. The random XORSAT model 

Let us now come back to our main topic, namely the behavior of random 
optimization problems, and consider them in the perspective of the mean field 
disordered models we have just discussed. We shall first emphasize the striking 
similarity between the energy function (2) of the XORSAT model and the one 
of the p-spin model (20): both are written as sums of products of Ising spin 
variables; for historical reasons the number of spins involved in each interaction 
is called k or p depending on the context. Apart from this minor conventional 
difference, the main discrepancy between the two cases is the structure of the 
interactions involved: in the random XORSAT problem there arc M = aN 
interactions with couplings of order 1, defining an Erdds-Rcnyi random hyper- 
graphs, while all the (^) possible couplings are present in (20), with individual 
strengths vanishing in the thermodynamic limit. Despite this difference both 
models are mean field, and share most of their phenomenology. In the XORSAT 
case there are two external parameters: the temperature T, and a "geometrical" 
parameter a, that controls the number of constraints put between the variables. 
It has been shown, in particular in [167, 168, 169], that the phase diagram in 
the (a, T) plane is divided in three regimes, separated by two transition lines 
otd(T) and a c (T), see Fig. 4 for a schematic representation. These two lines are 
the counterpart of the two transition temperatures Td and T c discussed above in 
the context of the fully connected p-spin model (that is recovered in the a — > oo 
limit). In the high temperature/ low a phase, the configuration space is well 
connected and ergodic; in the intermediate phase it becomes split in an expo- 
nential number of pure states, yet no singularity appears in the thermodynamic 
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Figure 4: A sketch of the phase diagram of the random XORSAT model in the (a, T) plane. 

functions on the line ad(T); the thermodynamic phase transition lies on the 
condensation line a c (T). 

We shall now concentrate on the zero temperature limit of the XORSAT 
model, i.e. on the properties of its ground state configurations, that are ob- 
viously the most relevant ones when the model is viewed as an optimization 
problem. The two transition lines have finite limits when T — > 0, that we shall 
denote ad and a c ; their expression as a function of k can be found in [167, 168]. 
It turns out, in this particular case, that a c = a s , where a s is the SAT-UNSAT 
threshold defined in Eq. (17). The dynamic transition ad is called clustering 
transition in this context, and similarly the pure states introduced above be- 
come clusters of solutions, i.e. sets of close-by solutions, well separated from 
each other. 

The XORSAT problem has some specific features that allow for an explicit 
definition of clusters which can be explained as follows [167, 168]. Consider an 
arbitrary XORSAT formula, and suppose that one of the variables <Xj appears 
in a single interaction, call it a. A moment of thought reveals that the formula 
is satisfiable if and only if the formula, with the interaction a removed, is sat- 
isfiable. One can iterate this process and reduce further the formula, removing 
at each step the interactions in which appears a variable of degree 1. At the 
end of this "leaf-removal" process one ends up with a reduced formula called 
the 2-core of the original one. Two cases can occur: either the 2-core is empty, 
and then the original formula is obviously satisfiable. One can assign satisfying 
values for the variables in the last removed interaction, and then reintroduce the 
interactions in the reverse order of the removal, using the fact that at each step 
at least one variable (the leaf) can be freely chosen to satisfy the re-introduced 
interaction. This case occurs with high probability when a < ad, and one can 
show [170, 171] that all the solutions that can be constructed from the free 
choices arc in some precise sense close one to each other. On the other hand, 
when a > ad the 2-core contains typically an extensive number of variables and 
interactions. This reduced formula, in which all variables are involved in at least 
two interactions, goes from satisfiable to unsatisfiable at the higher threshold 
a s [167, 168]. Let us consider the intermediate regime a € [ad,a s ], where the 
reduced formula on the 2-core is non-trivial but still has some solutions. A very 
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important point is that two distinct solutions a and a' of the 2-core formula are 
far away from each other, in the Hamming distance sense (the number of differ- 
ent variables between them). Indeed, because of the form of the constraints of 
the XORSAT problem, each interaction must contain an even number of spins 
i with Cj 7^ a[. In other words, as the 2-core does not contain leaves, a loop 
of disagreeing spins between a and </ has to be closed. As the random graphs 
are locally tree-like such a loop has necessarily a length diverging with N in 
the thermodynamic limit. Now, from every solution of the 2-core reduced for- 
mula one can construct different solutions of the full formula, by reintroducing 
the interactions in reverse order, as explained above. All these solutions that 
emerge from the same seed, i.e. from the same solution of the 2-core, will be said 
to belong to the same cluster (or pure state); then one realizes that solutions 
inside one cluster are close to each other, while solutions belonging to distinct 
clusters are necessarily separated by a large Hamming distance. To finish the 
connection with the phenomenology of the p-spin fully connected model, let us 
call e NS ( a " > the number of solutions of the 2-core formula for a random instance 
with parameter a; S(a) is defined in the interval [ad, ol c — a s ], precisely like the 
intermediate regime of temperature [T c ,Td] of Sec. 3.4. Moreover the complex- 
ity (or configurational entropy) S that counts the number of relevant clusters 
vanishes at the transition a c , similarly to the condensation on a sub-exponential 
number of pure states for T < T c . 

Let us summarize the main messages on the properties of random CSP that 
should be drawn from this particular case (see Fig. 5 for an illustration). For 
low values of the control parameter a the exponentially many solutions are 
spread in the whole configuration space, and close-by one to the other. In- 
creasing a there appears a clustering transition at ad, after which there are 
still exponentially many solutions, yet they are grouped in clusters of close-by 
configurations, the clusters being separated one from the other; in this regime 
the complexity or configurational entropy S counts the exponential number of 
clusters. The total entropy density of solutions, stot, is the sum of the com- 
plexity X and the internal entropy density s of each cluster, which is here the 
same for all clusters: Stot(a) = S(a) + s(a). For even larger values of a the 
satisfiability transition a s is due to the vanishing of S, i.e. the disappearance 
of the clusters of solutions; the last clusters that disappear can still contain an 
exponential number of solutions, i.e. the internal entropy density s(a) can be 
finite right at a s . The complexity £(a) and the total entropy Stot(a) are re- 
ported in Fig. 5 for 3- XORSAT on an Erdos-Renyi graph. The plot shows that 
indeed stot(«s) = s(a s ) is finite at a s for this model. This is important because 
it shows that typical instances have an exponential number of solutions even at 
the SAT-UNSAT transition, hence instances with a unique solution (that are 
particularly important for the analysis of the quantum adiabatic algorithm) are 
everywhere exponentially rare in this model. We will come back to this point 
in Sec. 3.9. 

We should also emphasize that XORSAT exhibits some specific features 
that are not shared by more complicated random CSP like fc-satisfiability or 
g-coloring. In particular all the clusters of XORSAT contain exactly the same 
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Figure 5: (Top panel) A sketch of the configuration space of random XORSAT problems. 
For low values of a the solutions, represented as black dots, are evenly spread on the N- 
dimensional hypercube. In the intermediate regime they are grouped in clusters, symbolized 
by the circles. For a > a B there are no more solutions. 

(Bottom panel) The total entropy of solutions stot(c) and the complexity (or entropy of 
clusters) for the 3-XORSAT problem on an Erdos-Renyi graph [167, 168]. The inset is 

a zoom of the region close to ay and a B . 
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number of solutions, because of the linear structure of the set of equations mod- 
ulo 2 it encodes. In general there are clusters of different sizes, and because of 
these fluctuations the condensation and satisfiability threshold do not coincide, 
as will be discussed further below. 

3.6. The random subcubes model 

The more complex phenomenology of random fc-SAT and g-COL has been 
first unveiled with rather intricate computations based on the cavity method, 
that we shall review in Sec. 5.1. For pedagogical reasons we shall first ex- 
plain this phenomenology using a toy model introduced in [172], the Random 
Subcubes Model (RSM), which is a non-local model (in the sense of Sec. 2.2.2) 
similar to the REM. The main new ingredient that is introduced in the RSM (to 
mimic fc-SAT and g-COL) is a distribution of clusters of different sizes. While 
in the XORSAT problem, for a fixed a, all clusters contain the same number of 
solutions, in the RSM it is assumed by construction that each cluster contains 
a different number of solutions, given by e Ns . Similarly, the number of clusters 
of internal entropy density s is given by e N ^ s ' . Hence the complexity is here a 
non-trivial function of s, like in the p-spin model discussed in Sec. 3.4, and not 
just a number as in the XORSAT problem. As discussed in Sec. 3.4, the fluctua- 
tions of internal entropy of clusters have an important consequence: they induce 
a new phase transition characterized by a condensation of the Gibbs measure 
onto a small number of clusters. Because in the RSM clusters are uncorrelated, 
all of its properties, and in particular the condensation transition, can be ex- 
tracted with much simpler computations than in random fc-SAT and g-COL. In 
addition, the RSM will be very useful in the quantum setting for understanding 
the effect of quantum fluctuations on random optimization problems. 

For all these reasons, we will discuss the RSM in more detail than we did 
for the previous models. In this section we explain the classical version of 
the random subcubes model [172], its quantum extension [47] being treated in 
Sec. 4.4. It will be useful for the discussion of Sec. 4.4 to define directly the 
model in a quantum notation, so we will do this here. 

3.6.1. Definition of the model 

The RSM distinguishes configurations that belong to a set of low energy 
clusters from those that belong to the remaining set of high energy configura- 
tions. It is defined as follows. Consider the Hilbert space H of N spins 1/2 
(qubits), in the basis of the Pauli matrices a?, \g) = \a±, • • • , o~n)- A cluster A 
is a subset (subcube) of the Hilbert space 

A = {\a)\Vi : aiCirf} , (27) 
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Figure 6: Pictorial representation of the different phase transitions in the set of solutions of 
the random subcubes model [172]. 



where irf are independent random sets defined as follows: 

{— f with probability — 1 

p > <Ji is "frozen" in cluster A, 

1 with probability — 

{1,-1} with probability 1—p <7j is "free" in cluster A. 

(28) 

Thus, with probability p the variable i is frozen in A and with probability 1 — p it 
is free. With this definition the number of states, i.e. classical configurations, in 
a cluster A is a random variable equal to 2 Nsl -- A \ where Ns(A) is the number of 
free variables and we call s(A) the internal entropy of a cluster (for convenience 
in this section we use log 2 to define entropies). We next define a set S as the 
union of 2 w ( 1_a ) random clusters, and its total entropy stot' 

2 JV(l-cv) 

S = |J At stot = ^ log 2 \S\ . (29) 

i=l 

The parameter a here is analogous to the density of constraints in CSP. The 
probability p that a variable is frozen instead plays the role of the clause size k 
in fc-SAT or the number of colors q in the g-coloring problem. 

For each cluster A we assign a Hamiltonian H A = Neo(A) ^2 aeA \s)(s.\ with 

e o(^4) > and a "penalty" Hamiltonian Hy = NVj^afs \s)is.\ which describes 
the classical energy of states not belonging to S. The problem Hamiltonian 
Hp = Hy + a Ha is of course diagonal in the basis |cr), and the associated 
cost function E(a) = (a\Hp\a) is equal to NJ2 A . a(£A eo(A) if a <G S and NV 
if a $l S. With these definitions we wish to interpret the states in S as "local 
minima" of Hp and the others as "excited states" . A sharp distinction between 
them can be obtained by sending the positive constant V to infinity, a choice 
that we adopt in this section. In Sec. 4.4 we will consider also a finite V, but 
we will always assume that V ^> max^ eo(A). 

3.6.2. Clustering 

To begin the discussion, we shall briefly characterize the structural changes 
in the set S when a is varied, which are shown in Fig. 6 and have been derived 
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in [172]. We will assume here for simplicity that all the clusters have eo(A) = 0. 
The set S is then the set of the "solutions" of the problem (or of the ground 
states of zero energy). We will show that the transitions ad, a c and a s outlined 
in Sec. 3.1 can be very precisely defined in the RSM. 

We will extensively make use of two well known results of probability theory 
called the union bound and the Chebychev inequality. In fact, the properties of 
iS can be traced back to probability statements concerning 2 Nb events £j, each 
having probability "P(£,:) = 2~ Na , for some a and b. Under these conditions the 



union 



bound states: 



2 Nb 

<^V{£i) = 2 N{ - b - a \ (30) 

i=l 



which implies that when a > b the probability V(Ui£i) is exponentially sup- 
pressed in the size of the system N. When the events are independent, the 
number of true events Af is a random variable with a binomial distribution with 
(TV) = 2 Nl - b '^ and (Af 2 ) = 2 Ni ~ b - a \l - 2~ Na ). Then for arbitrary small e one 
can apply the Chebychev inequality: 

i^-WLy m < 1 (3D 



(AT) ) ~ (AO 2 

which ensures that when a < b, Af is self-averaging in the large N limit, i.e. the 
average is exponentially large and concentration around the average Af ~ (Af) 
is found. 

As a first application, we consider the number of clusters Af(s) of entropy s. 
Because frozen variables are chosen independently, we have 

V(a(A) = 8)=^p»V-\l-p) If '. (32) 

Hence Af(s) follows a binomial distribution with parameter V{s) and 2 Ar ' 1 ~ Q ) 
terms, and its average, at the leading exponential order, is (Af(s)) = 2 Ar ( 1 ~ a - ) 7- > (s) 
2 NS ( S ) with the complexity S(s) defined by 

Z(s) = l-a-D(s\\l-p) , 

D(x\\y) =xlog 2 (x/y) + (1 - z)log 2 [(l - x)/{\ - y)} . 

In the region s € (s m i„, s max ) where S(s) > 0, we can apply the Chebychev 
inequality to show that Af(s) concentrates around the average when N — > oo. 
In the region where S(s) < 0, we can apply the union bound to show that with 
probability 1 there are no clusters of entropy s when A^ — > oo. 

Next, we apply similar arguments to identify the following changes of the 
structure of the space of solutions when a is varied (see Fig. 6). We will only 
give brief sketches of the proofs; the reader is referred to [172] for further details. 

• For a < «d = l°g2(2 — p), each state \a) belongs to an exponential number 
of clusters and S = H. For a > Qd a random state does not belong to S 
with probability 1 when A^ — > oo, thus S ^ % and stot < 1- 
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Proof: The probability that a configuration \a) belongs to a cluster A is 
V(\q) S A) = (1 - \ ) N and 



2 

P 



N 2 «(l- Q ) 

V(\a) i S) = [1 (l |j . (34) 

Then from the union bound if a < ad — log 2 (2 — p): 

v(s + u) = v(u^\s_) i s) < 2 N e -2^-^ ^ (35) 

which implies that all states are in S and Stot = 1- For a > ay, from 
Eq. (34) we get V(\a) (£ S) 1. Thus S ^ H and s to t < 1- 

For a > a scp — 1 + log 2 (l —p 2 /2)/2, the clusters are well separated, in the 
sense that with probability 1 for N — > oo the Hamming distance (minimal 
number of different spins) between any two clusters is of order N. 

Proof: Wc note that V(A n A' ^ 0) = (1 - ^-) N . Then we can apply the 
union bound over all possible intersections in the set S 

i / 2 \ N 

V(U ll (A. l nA 1 ^H}))<-2 N ^- a \2 N ^ ^0 (36) 

for a > a scp . This means that with probability 1 when N — > oo the clus- 
ters are disjoint, i.e. their Hamming distance is strictly positive. The prob- 
ability to hnd clusters at distance x is finite only when x = 0{N) [172]. 

For Qd < a < a c = p/(2 —p) + log 2 (2 —p) most of the solutions belong to 
one of the exponentially many clusters of size s*, with S(s*) > and s* <G 
(•Smin, s mal ). On the contrary when a > a Cl s* = s max and most of the 
solutions belong to the largest clusters whose number is sub-exponential 
in N because E(s max ) = 0. 

Proof: One can compute the total number of states in S by observing that 

\S\ = 2 Wstot ~ Yl 2NS{A) ~ / SmOX ds 2 m{s)+s] , (37) 

therefore s tot = max se [ Smm , Smax ] P(s) + s]. Studying the function S(.s) + s 
it turns out that up to a c its maximum value, dominating the saddle point 
in the integral, is taken inside the allowed interval and thus S(s*) > 0. 
When a > a c instead the maximum is achieved at the boundary of the 
interval, implying W(s*) = 0(1). 

Finally, for a > a s = 1 there are no more solutions. 

Proof: This follows trivially from the definition of the number of clusters, 
equal to 2 N( - 1 -°\ Then for a > 1 there are no more clusters and the 
set S is empty. In the language of random CSP, a s corresponds to the 
SAT-UNSAT transition. 
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Note that in this particular model the entropy has a singularity at ad , which 
is not present in local random CSP. From the dynamical point of view what char- 
acterizes ad is that for a > ad there is "ergodicity breaking" in the sense that a 
local random walk over solutions starting in one cluster takes an exponentially 
long time to reach another cluster [172]. 

3.6.3. The partition function at finite temperature 

Similar results can be obtained when the clusters have a distribution of 
energies [172]. Let us assign to each cluster an i.i.d. random energy eo € [0,e m ] 
in such a way that the total number of clusters of energy eo is 2 N ^ 1 ~°^ 9 ^ e °\ 
with (?(eo) an arbitrary increasing function of eo, because it is reasonable to 
assume that the number of clusters increases with energy. Then, the above 
arguments can be easily generalized for each level of energy eo. Following the 
same reasoning that leads to Eq. (33), the number of clusters of energy eo and 
entropy s is 2 N ^ (e °' s \ with 

E(e ,s) = (l-a)p(eo)- D(s\\l -p) , (38) 

and is positive in an interval s £ (s m j n (eo), s max (eo)). Of particular interest is 
the computation of the partition function at finite temperature, that replaces 
Eq. (37) and reads 

/•e m /■Smax(eo) 

z = Y^ 2 Ns{A) e -PNe {A) ^ / ^1 dg 2 N[S(e ,s)+ S -^ log 2 e] ^ (3g) 

A JO A min (e ) 

and that can be evaluated by a saddle point. When the saddle point values of 
eo and s reach the boundary of the integration interval a condensation transi- 
tion happens, on a line a c (T). Moreover, at each level of energy, the previous 
analysis of the structure of the union of clusters can be repeated, and the same 
transitions happen at energy- dependent values ad(eo), a scp (eo). Because for 
each temperature a unique value of eo dominates the partition function, these 
can be converted in lines ad(T), a sep (T). It is easy to show that all transition 
points increase with T, like in the XORSAT case (Fig. 5). 

3. 7. The space of solutions of random constraint satisfaction problems 

By means of non-local toy models such as the REM (Sec. 3.3) and RSM 
(Sec. 3.6) we built a lot of intuition about the different transitions that happen 
in spin glass models. The analysis of the XORSAT problem (Sec. 3.5), that 
can be carried out in a relatively straightforward way using rigorous methods, 
also illustrated the emergence of clustering in the simplest random CSP. Hav- 
ing now understood the kind of transitions to be expected, one would like to 
make precise computations in other local random CSP such as /c-SAT or q-COL. 
Unfortunately, this turns out to be a more difficult task and requires the intro- 
duction of sophisticated statistical mechanics methods, in particular the cavity 
method [137]. These will be introduced in Sec. 5. 
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Figure 7: Sketch of the space of solutions — colored points in this representation — in the 
g-coloring problem on random graphs when the connectivity c is increased [130, 173]. (i) At 
low c, all solutions belong to a single cluster, (ii) For larger c, other clusters of solutions 
appear but a giant cluster still contains almost all solutions, (iii) At the clustering transition 
Cd, it splits into an exponentially large number of clusters, (iv) At the condensation transition 
c c , most colorings are found in the few largest of them, (v) The rigidity transition c r (c r < c c 
and c r > c c are both possible depending on q) arises when typical solutions belong to clusters 
with frozen variables (that are allowed only one color in the cluster), (vi) No proper coloring 
exists beyond the COL/UNCOL threshold c s . 



Before turning to quantum versions of spin glass models, we will give here a 
summary of all possible transitions that have been found in local random CSP at 
zero temperature, taking as an illustrative example the random g-COL problem 
with q > 4 colors (the q = 3 case being a bit particular) and a large Erdos- 
Renyi random graph whose average connectivity c = 2a increases. Different 
phases are encountered that we will now describe in order of appearance. The 
corresponding phase diagram is depicted in Fig. 7 [130, 173]. 

(i) A unique cluster exists: For low enough connectivities, all the proper 
colorings are found in a single cluster, where it is easy to "move" from one 
solution to another: for any given pair of solutions, one can construct a 
path of solutions that connects them, such that at any step along the path 
only a sub-extensive number of colors are changed. The total entropy of 
solutions can be computed and reads in the large graph size ./V limit: 

stot = log<?+ | log (l - -J . (40) 



(ii) Some (irrelevant) clusters appear: As the connectivity is increased, 
the phase space of solutions decomposes into a large (exponential) num- 
ber of different clusters. It is tempting to identify that as the clustering 
transition, but it happens that all (but one) of these clusters contain rel- 
atively very few solutions — as compared to whole set — and that almost 
all proper colorings still belong to one single giant cluster. Clearly, this is 
not a proper clustering phenomenon and in fact, for all practical purposes, 
there is still only one single cluster. Eq. (40) still gives the correct number 
of colorings at this stage. 
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(iii) The clustered phase: For larger connectivities, the large single cluster 
also decomposes into an exponential number of smaller ones: this now 
defines the genuine clustering threshold Cd- Beyond this threshold, a local 
algorithm that tries to move in the space of solutions will remain prisoner 
of a cluster of solutions. Interestingly, it can be shown that the total 
number of solutions is still given by Eq. (40) in this phase. This is because 
the free energy has no singularity at the clustering transition (which is 
therefore not a true transition in the sense of Ehrenfest, but rather a 
geometrical transition in the space of solutions). 

(iv) The condensed phase: As the connectivity is further increased, a new 
sharp phase transition arises at the condensation threshold c c where most 
of the solutions arc found in a finite number of largest clusters. From this 
point, Eq. (40) is no longer valid, because this is a genuine phase transition. 
The entropy is therefore non-analytic at c c and Eq. (40) becomes just an 
upper bound. 

(v) The rigid phase: Two different types of cluster exist: in the first type, 
that we shall call the unfrozen ones, all spins can take at least two different 
colors. In the second type however, a finite fraction of spins are allowed 
only one color within the cluster and are thus "frozen" into this color. It 
follows that a transition exists, that we call rigidity, when frozen variables 
appear inside the dominant clusters (those that contain most colorings). 
If one takes a proper coloring at random above c r , it will belong to a 
cluster where a finite fraction of variables is frozen into the same color. 
Depending on the value of q, this transition may arise before or after the 
condensation transition (a list of values can be found in [174, 173]). 

(vi) The UNCOL phase: Eventually, the connectivity c s is reached beyond 
which no more solutions exist. The ground state energy (sketched in 
Fig. 8) is zero for c < c s and then grows continuously for c > c s . The 
values c s computed within the cavity formalism are in perfect agreement 
with the rigorous bounds [120, 121, 122, 138] derived using probabilistic 
methods and are widely believed to be exact, although this remains to be 
rigorously proven (see [131, 132] for a proof that they are at least rigorous 
upper bounds). 

Notice that in specific models some of these transitions coincide. We have 
already seen in Sec. 3.5 that in XORSAT c r = Cd and c c = c s , therefore some 
of the phases above do not exist: all clusters are frozen, and the condensed 
phase does not exist. Another example is the 3-COL problem, which is peculiar 
because Cd = c c so that the clustered phase is always condensed. In view of this 
rich and model-dependent phase diagram, it is important to get an intuition on 
the meaning and the properties of these different phases. 

At this point, there are many questions one could ask. First of all: are 
these problems hard only close to the SAT-UNSAT threshold c s ? The answer is 
no: for instance in q coloring, when q is large, problems are easy (in this case, 
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the complexity is linear in the number of nodes) for almost every algorithm as 
long as c < Cd ~ q^ogq (to leading order) but suddenly very hard (so that 
no algorithm is known that performs provably in sub-exponential time in the 
number of nodes) if c > c<j. It is known, however, that there exist solutions 
up to c = c s ~ 2q log q. A similar problem appears in random fc-SAT between 
2 k \ogk/k and 2 fc log2 [175]. One could then conclude that the clustering above 
Cd is responsible for the hardness of the problem. Yet, for small enough q (e.g. 
q = 4), many algorithms are able to find solutions in the clustered phase at c 
much larger than Cd [126, 176, 177]. Why then are some problems hard and 
some easy? Does something else explain the sudden onset of hardness? 

Unfortunately, the answer to these questions is in large part still open. Yet, 
many interesting results on the connection between the above picture and algo- 
rithmic hardness have been obtained. For reasons of space, in the rest of this 
section we will focus in particular to simulated (thermal) annealing. 

3. 8. Efficiency of the simulated annealing 

It turns out that close to the satisfiability threshold c s finding the solu- 
tions to the problem becomes particularly hard and most algorithms suffer of a 
dramatical slowing down. Quite generally this phenomenon is attributed to the 
presence of many minima in the energy landscape and to the organization of the 
solutions in phase space. Simulated (thermal) annealing [86] is one of most fa- 
mous algorithms designed to tackle complex energy landscapes. Despite the fact 
that it only partially accomplishes this task as it actually fails when too many 
clusters dominate the partition function, it represented a true breakthrough 
in the domain and it is still exploited in many applications. The prescription 
of simulated annealing (Sec. 2.3.1) is to initialize the algorithm with a random, 
high temperature, configuration. Then, lower the temperature, eventually down 
to zero, in discrete steps according to an assigned protocol, and at each step, 
perform a given number of local movements in phase space -Monte Carlo stcps- 
in order to equilibrate at that temperature and use the last generated config- 
uration to initialize the search at the new temperature. Technically this is the 
implementation of a time dependent Markov chain. An implementation of sim- 
ulated annealing in continuous time is also possible. In this section we want to 
discuss in more details the relation between the structural transitions discussed 
in Sec. 3.7 and the performances of simulated annealing. 

3.8.1. Effects of the clustering transition on thermal annealing 

In order to discuss better the properties of thermal annealing, we need to 
introduce a temperature into the problem, as discussed in Sec. 2.3.1, and investi- 
gate the finite temperature phase diagram. The latter is sketched in Fig. 8 [178] 
for 5-COL with q > 4 on Erdos-Rcnyi random graphs as a function of average 
connectivity c. At high temperature the system behaves as a paramagnet in 
the language of magnetic systems. The clustering and condensation transitions 
extend in lines T d (c) and T c (c). On the contrary, the rigidity and SAT-UNSAT 
transitions exist only at zero temperature, because at finite temperature the 
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Figure 8: Typical phase diagram for a spin glass problem on a random graph. At Tj, local 
Monte Carlo dynamics becomes inefficient ("dynamical" transition [157, 169, 179]). At T c the 
system undergoes an equilibrium glass transition [180, 178]. eQs represents the ground state 
energy and is positive when the problem is not satisfiable. At zero temperature, the phase 
transitions of Sec. 3.7 are recovered. 

notions of "solution" and "frozen variable" cannot be defined: constraints can 
always be violated with some finite probability. 

The idea behind a simulated annealing in temperature is that thermally 
activated processes allow to overcome the energy barriers and at the latest 
stages the zero temperature dynamics converges towards the solution. However 
this is only true if the annealing is slow enough. Let us first try to apply, as 
in the quantum case, an adiabatic strategy. There exist rigorous bounds on the 
time needed for a thermal annealing to stay adiabatic, but they yield, in the 
worst case, an exponential time [181]. Indeed, on random CSP for c > Cd, it 
turns out that equilibration is an exponentially hard task when the clustering 
temperature 'Td(c) is reached. Below Td(c) the dynamics falls out-of-equilibrium 
[179] unless one is ready to wait for exponentially long times. This rigorous 
result can be intuitively understood by the following three arguments. First, 
the probability to overcome the barriers between the pure states of the Gibbs 
measure is exponentially small in the size of the system, hence these states trap 
the dynamics for an exponentially large time. Secondly, even if one can go out 
of a pure state, there are exponentially many of them, and it thus takes again 
some time to find the equilibrium ones [182, 183]. To further complicate the 
problem, a third effect exists: even if one manages to equilibrate the system 
at a given temperature T, this is not really useful because the pure states that 
dominate the partition sum at any T' < T arc completely different ones, so 
that the hard equilibration work has to be entirely redone from scratch as soon 
as the temperature is slightly changed. This is an effect called temperature 
chaos [184, 185, 186, 187, 188] which is the classical analog of the quantum level 
crossing that will be discussed in Sec. 4.4. The combination of all these effects 
is behind the exponential hardness of an adiabatic thermal cooling. 
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3.8.2. Non-adiabatic thermal annealing 

If one, however, is not interested in being adiabatic and instead is just inter- 
ested in the final configuration reached, then the situation is different: it may 
be possible to reach a zero energy state at the end of the protocol while being 
out of equilibrium at intermediate stages of the annealing. This is observed for 
instance in random coloring [189]. There are also many random walk algorithms 
that are similar to Monte Carlo Markov chains but do not satisfy the detailed 
balance condition: these can find a solution of random CSP in some part of the 
clustered phase [190, 177]. This seems at first contradictory with the previous 
section. However, the problem of finding a solution is different from sampling 
solutions uniformly, which is what is achieved by an adiabatic cooling. If one 
just wants to find a solution, it is often possible to succeed up to much larger 
connectivities than the clustering one. In fact, it seems from empirical evidences 
(see the discussion in [187, 188, 177]) that the moment the problems become 
truly hard is when the rigidity transition is reached, or more precisely when all 
solutions belong to frozen clusters. 

This can be understood in terms of energy landscape: below the clustering 
temperature, one is trapped into a single pure state and cannot visit the whole 
space. However, the energy of the pure state one is trapped in can be lowered 
when the temperature is reduced, and maybe even go to zero when T — >• 0. As 
shown in [187, 188, 177] such "canyon-like" states, that reach the zero energy 
configurations, become rare after the rigidity transition. In this case, one needs 
again to visit many states until a good one is found, and one is back to the 
situation discussed in the previous section. In summary, if one is ready to forget 
about adiabaticity, annealing and other strategies can be applied and work well 
up to connectivities larger than cj, but fail when rigidity is met. According 
to this empirical analysis, the clustered and frozen region of the phase diagram 
contains the hardest possible instances of random CSP. 

3.8.3. Existence of good paths for classical annealing 

We have discussed so far annealings in temperature, but other annealing 
protocols can be used: for instance, an annealing starting from a large magnetic 
field and reducing it to zero, or any path in a temperature-field phase diagram 
that ends at zero field and zero temperature. Of course, the fact that a thermal 
annealing is inefficient does not imply that all possible paths are, eventually, 
inefficient. Indeed, it is easy to see that there are paths that will make a classical 
annealing work in a polynomial time, however hard the problem is. 

Consider for instance the following setting: we take an instance of a CSP 
defined on Ising spin variables, and one of its solutions r = {n}, where n is 
the value of variable i in the solution. Consider now a protocol in which the 
following term is added to the cost function: 



The "local magnetic field" hi = hri points each spin in the direction of the 
correct value in the chosen solution. We can now perform a simulated annealing 




(41) 
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by starting at large T with finite h and then reducing progressively T and h 
to zero. Upon cooling in this field the final configuration will be the solution 
associated with the field. One might (rightfully) argue that using a protocol 
that knows the solution of the problem is equivalent to cheating. Indeed the 
fact that such a protocol exists is not very useful, because finding the right field 
is equivalent to finding a solution. The problem is then: how difficult is it to 
find a "good" annealing path? 

The message here is that the mere fact that some efficient annealing protocol 
exists is not conclusive. More generically, if one studies classical algorithms to 
solve CSP, there is always a good one (for instance, if one initialize a random 
walk algorithm in a solution) but it is, of course, hard to find it in general. To 
prove rigorously this statement is however difficult; this is basically proving that 
P is not NP. In short: the question is not to decide whether there is a classical 
protocol able to find a solution quickly: for a given problem, such protocols 
always exist. The question is rather how to find them. The only thing one can 
do is to consider a given protocol and study it: to the best of our knowledge, 
the random instances in the clustered and entirely frozen region are extremely 
hard — and can be considered as some of the hardest representatives of the NP 
class — and there has been no practical way to solve them generically. 

One will have to keep these considerations in mind in the analysis of the 
quantum annealing. Proving that efficient quantum paths exists is not enough 
from an algorithmic perspective: one really needs to be able to construct these 
paths explicitly. After all, we are interested here in constructing a specific 
algorithm (a specific annealing protocol) and prove its efficiency (or inefficiency). 

3.9. Generating USA instances (locked CSP) 

Instances with a Uniquely Satisfying Assignment (USA) are very practical 
in quantum studies, as the absence of degeneracy of their ground state allows to 
unambiguously define their gap. Generating a problem that has only one single 
solution is, however, not so easy a priori. Actually, even certifying that an 
instance is USA requires to count all of its solutions, which is an exponentially 
hard task (therefore the decision problem of whether an instance is USA or not 
is in general not in NP). 

Moreover, the usual protocol, used since [20], consists in generating many 
instances of a given problem (for a fixed size) , finding all the solutions of each in- 
stance using a complete solver, and then keeping only the instances with a single 
solution. It turns out that the latter is a very difficult task, that is possible only 
for small sizes. In fact, in most random CSP, there are generically exponentially 
many solutions, so that one needs to generate exponentially many instances to 
find the rare ones with a USA and in the end generating the instances in this 
way is at least as hard than finding their solution. Fortunately, there is a way 
around this problem, using the so-called "locked problems" at the satisfiability 
threshold, or using a hidden assignment in their UNSAT phase. 
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3.9.1. Locked problems at the SAT- UN SAT threshold 

The concept of "locked problems" was introduced in [191]. This is a very 
broad class of CSP where i) if one changes one variable (and only one) in a 
given clause in a SAT configuration, then the system becomes UNSAT and 
ii) each variable is involved in at least two clauses. If one wants to keep the 
configuration SAT, then one has to flip a variable in a neighboring clause, and 
so on, until a loop is found and the chain can be closed. This means that, on 
a random graph that is locally tree-like, in order to go from one solution (or 
satisfying assignment) to another one, it is necessary to flip at least a closed 
loop of variables in the factor graph representation, which typically involves 
O(logiV) changes. 

There are many such models, such as the XORSAT, l-in-3 and l-in-4 SAT 
problems (Sec. 2.1.1) defined on the ensemble of regular random graphs [191]. 
Generically, the phase transition in locked models are the same as in the XOR- 
SAT problem on an Erdos-Rcnyi graph, as illustrated in Fig. 5, but with some 
important differences: 

• When a is small enough, the set of solution is "nearly" connected: by 
flipping O(logiV) variables, one can visit all possible solutions starting 
from a particular one. In fact, if one allows also to visit "quasi solutions" 
with O(l) cost, single flip spins are enough to visit the whole space of 
solution 1 . We thus say that the space of solutions is made of a single, 
unique cluster. 

• When a > ad, the set of solution undergoes a clustering transition and 
splits into an exponential number of components, separated by O(N) flips. 
This is the clustering transition. The crucial characteristic of locked model 
is that each cluster contains a single solution (the circles in Fig. 5 contain 
only one black dot in this case). Therefore, the internal entropy of clusters 
vanishes, s(a) = 0, and the total entropy coincides with the complexity, 
stot(a) = S(a). 

• For a > a s , there are no solutions anymore. As usual, the complexity van- 
ishes continuously at the SAT-UNSAT transition a s . Because s to t(aO = 
S(a), the total entropy also vanishes at the transition: s to t(o; s ) = 0. 

There are two consequences of these properties that are important for the 
present discussion. First, these models are hard to solve not only at the thresh- 
old, but also in the full clustered phase (see [193, 191]): here, the clustering and 
rigidity transitions coincide and the entire clustered phase is always associated 
with hard instances. In these models we thus have a clear and well defined link 
between the presence of a phase transition and the computational complexity. 
The same conclusion is valid in particular for the XORSAT problem on random 
regular graphs. This problem is in P (while locked problems are generically 



1 For expert readers: this is the difference between reconstruction and small noise recon- 
struction [191, 192]. 
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NP-complete) , as it can be solved by Gaussian elimination. If one, however, 
decides to forget about this information, it is a very difficult problem. In fact, 
the hardest instances of SAT problems to this day are obtained by generating a 
XORSAT formula and by adding to it a bit of non linearity (such that Gaussian 
elimination cannot be used, see [194, 195]). 

There is another reason why these models are useful in quantum annealing 
studies. When working exactly at the satisfiability thresholds a B , the entropy is 
exactly zero as we discussed above: hence, USA instances appear with a finite 
probability, making them very easy to generate (see Sec. 6.2.3 and [48]). 

3.9.2. Planted locked models 

There exists another way to generate USA instances, not only with finite 
probability, but with a probability going to one in the large size limit: planting 
a solution in the UNSAT phase of locked models. This was proposed in [196] 
and studied with statistical physics methods and rigorous mathematical proofs. 

The idea is to create both an instance and a solution by first assigning a 
configuration to all variables, and then choosing only constraints compatible 
with this configuration. This creates instances from the so-called "planted" 
ensemble. As shown in [196], for random locked CSP models, instances from the 
planted ensemble have with high probability a single satisfying assignment (or a 
pair of them if a global symmetry is present) beyond the satisfiability threshold. 
This allows to create USA instances of any size at zero computational cost. 

The question is whether the instances created in this way are difficult. This, 
again, depends on the average density of constraints [196]: an easy-hard-easy 
pattern for finding a solution appears in the planted ensemble as the constraint 
density is increased. The boundaries of the hard phase are given by the clus- 
tering transition on one side, and by another transition (called the threshold 
for the robust reconstruction [196]) on the other side. Between these two tran- 
sitions, hard instances with USA are generated. Again, XORSAT on random 
regular graphs is particular in the sense that planted instances remains hard for 
arbitrary large connectivity 

We shall therefore discuss in a lot of details the quantum XORSAT model 
in Sec. 6, using its double status as a hard benchmark and as a simple model 
to generate USA instances. 

4. The low energy spectrum of quantum spin glasses 

As described in Sec. 2 the Quantum Adiabatic Algorithm (QAA), or quan- 
tum annealing, is a procedure designed to find the ground state of a classical 
Hamiltonian. From the point of view of the computational complexity theory, 
random optimization problems are thus natural benchmarks for this algorithm. 
We described in Sec. 3 the rich phenomenology of these classical disordered mod- 
els, and explained how the study of mean field spin glasses allowed to understand 
them. We shall now progressively turn towards the study of quantum optimiza- 
tion problems, i.e. classical disordered models with some non-commuting term 
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(for instance a transverse field) inducing quantum fluctuations. As in Sec. 3 we 
will start by discussing simpler models for pedagogical reasons. 

The quantitative assessment of the performances of the QAA, as discussed 
in Sec. 2.3.2, requires a detailed understanding of the low energy spectrum of 
disordered quantum Hamiltonians. Its time complexity is indeed directly related 
to the square of the inverse of the minimum gap A m ; n of the Hamiltonian along 
the annealing path. On the basis of complexity theory we are then particularly 
interested in discriminating between Hamiltonians whose minimum gap vanishes 
polynomially and those for which A m ; n is expected to be exponentially small 
in N. 

It is well known that the gap of the Hamiltonian vanishes upon increasing N 
in correspondence with a quantum phase transition [33, 38, 48, 39, 197]. How- 
ever, the reverse is not true: the gap might vanish even without an underlying 
quantum phase transition. Indeed, some model Hamiltonians display phases 
such that the gap is everywhere exponentially small in N, due to a continuum 
of level crossings [47] . Interestingly, this was associated to a kind of Anderson 
localization phenomenon in phase space [40]. In these cases, the vanishing of 
the gap in the thermodynamic limit is not associated to a singularity in the 
ground state energy [47]. 

In this section, we will present a review of quantum phase transitions in 
several simple Hamiltonians, and discuss the corresponding scaling of the gap. 
We will discuss how level crossings can be induced by disorder, and how their 
accumulation can result in a complex spin glass phase where the gap is every- 
where exponentially small. We will keep the discussion informal, and focus on 
toy models. A more precise discussion on realistic optimization problems will be 
presented after the methods to study such complex phenomena will have been 
introduced in Sec. 5. 

For concreteness in this section we shall only consider models of quantum 
spins 1/2, with Hamiltonians H made of two terms. The first is diagonal in 
the eigenbasis of the af operators, it encodes the problem to be solved and 
we will refer to it as Hp. The second term induces quantum fluctuations of 
strength T and we will call it THq, in such a way that the total Hamiltonian 
is H = Hp + THq. For concreteness, we will consider as a quantum term a 
transverse field, Hq = — 5Z i=1 Z?f . Hence we denote here V the strength of 
this transverse field, and when speaking of an annealing it is understood to be 
from r = oo down to V = 0. The connection with the notations of Sec. 2 is 
easily made: the problem Hamiltonian Hp corresponds to H[ , the quantum Hq 
corresponds to H^, and the correspondence between the interpolation parameter 
s and r is r = (1 — s)/s. Quantum fluctuations different from a transverse field 
are expected to have similar effects, provided they are simple enough and in 
particular they do not contain detailed information on the classical part of the 
Hamiltonian (for instance, a hopping quantum term —t^2uj) {&f&j was 
considered in [198] and led to similar results). 
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4-1- Second order transitions 
4-1.1. Ordered models 

A well known example of a system without disorder that exhibits a second 
order phase transition associated to a polynomially vanishing gap is the one 
dimensional Ising ferromagnetic chain in a transverse field: 

N N 
i=l i=l 

This Hamiltonian is integrable and can be solved completely [33] . Generically, in 
finite dimensional ordered models, the gap of the system, in the thermodynamic 
limit, vanishes when one gets close to the quantum phase transition as t zv , 
where e is a measure of the distance to the transition, v is the exponent for the 
divergence of the correlation length and z the dynamic exponent. Right at the 
transition the gap vanishes only in the thermodynamic limit, it thus scales with 
the size N of the system, for instance as 1/N in the one dimensional case. As 
our goal is to study mean field models, for the motivations explained in Sec. 3, 
here we want to mention also the Curie- Weiss mean field Hamiltonian [199, 200, 
201, 202]: 

N N 

g = -2ivE^r r E^' ( 43 ) 

z,j = l i—1 

This Hamiltonian is easily solved by the mean field construction, that amounts 
to replace one of the a z by its average m; one then obtains a single site Hamil- 
tonian H = —Jma z — Ta x and the magnetization is computed self-consistently, 
leading to the mean field equation 

to = Jm tanh((3y/j 2 m 2 + T 2 ) . (44) 

v j 2 to 2 + r 2 

At zero temperature, this leads to a phase transition between a paramagnetic 
(to, = 0) phase at T > J and a ferromagnetic phase with magnetization to = 
y/l — (r/J) 2 for r < J. The order parameter m is continuous at the transition 
and the ground state energy as a function of V has a singularity in the second 
derivative. Hence the transition is of second order. It is possible to show that 
the gap of the Hamiltonian vanishes polynomially (more precisely as iV -1 / 3 ) 
at the phase transition point r = J [199, 200, 201, 103], therefore a quantum 
annealing can find its ground state in polynomial time. 

4-1-2. Disordered models 

Let us now consider what happens when disorder is introduced in the Hamil- 
tonians described above. In the unidimensional case, Eq. (42) now becomes: 

N N 
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Consider for instance the case where the Ji and I\ are i.i.d. random vari- 
ables, uniformly distributed in [0, J] and [0, V] respectively (negative couplings 
or transverse fields can be eliminated through a simple redefinition of the spins). 
Note that in the classical limit T = the ground state is <Tj = 1 or <x; = — 1, 
hence there is no frustration in the model. 

This random model has been extensively studied by means of the renor- 
malization group by Fisher [203]. For our purposes, the main results are the 
following: 

• A quantum critical point is present at T = J, such that for T < J the 
system is ferromagnetic while for T > J it is paramagnetic [203] . 

• At the transition point T = J, the typical gap is exponentially small: 
A = e~~ 9 ^ ', where g has a finite probability distribution over disorder 
(hence the distribution of A is very broad) [204, 205]. 

• Above the transition, V > J, the gap is typically of order one [204, 205]. 

• Below the transition, T < J, the gap is exponentially small, but its 
logarithm is concentrated around its average: it is distributed accord- 

2 

ing to a Gaussian distribution with log A cx N, and (log A) 2 — log A oc 
N [204, 205]. 

• The exponentially small gaps for T < J are due to crossings within differ- 
ent low energy levels, due to the fact that these levels are localized (in a 
sense that will be made more precise below) [102, 39, 40]. 

Hence, a quantum annealing will typically encounter a gap A ~ exp(— g^/~N) at 
the transition T = J, followed by a series of gaps A ~ exp(-g'N) for V < J, and 
will typically require an exponential time to find the ground state. However, 
it has been shown that the residual energy per spin after a quantum annealing 
over a finite time r (as well as after a simulated annealing) goes to zero when 
t — y oo [102, 206]. This implies that although finding the ground state is 
exponentially hard in N, finding a state whose energy per spin coincides with 
the one of the ground state is indeed quite easy (it can be done in polynomial 
time) 2 . 

The mean field model that corresponds to Eq. (45) is the quantum Shcrrington- 
Kirkpatrick (SK) model, which can be obtained either from the quantum Curie- 
Weiss model of Eq. (43) by including randomness in the interaction couplings, or 
from the classical SK model defined in Eq. (19) by the addition of a transverse 



2 It was shown in [102, 206] that the total residual energy goes as AE ~ 7V/(logr)« at large 
times. This implies that finding the ground state energy (i.e. finding a state with energy AE 

of order 1) requires a time r ~ e N ^ . However, if one is only interested in finding the ground 
state energy per spin, it is enough to require that AE grows slower than N. For instance one 
can choose AE = N/ log N and in this case a time r ~ N 1 ^ is enough. 
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field: 

N 

H=-Y / J ij a^-rY i df . (46) 

Here, the Jy are quenched i.i.d. Gaussian variables, with zero average and 
variance Jf- = J 2 /N; hence, the couplings can be positive or negative, leading 
to frustration. This model has a quantum critical point at T = J, separating 
a spin glass phase at T < J from a paramagnetic one at Y > J [207]. The 
analysis of the gap in this model is however more difficult, and it has not yet 
been performed to our knowledge. However, it has been shown [208] that the 
whole spin glass phase at T < J is gapless for N — >• oo. Because the spin glass 
phase is characterized, in the classical limit r = 0, by many almost degenerate 
low energy minima, it is very natural to expect, like in the finite dimensional 
case, the existence of level crossings between the energy levels corresponding to 
these minima when quantum fluctuations are switched on [102, 39, 40]. These 
should lead to an exponentially small gap in the whole spin glass phase, as we 
will discuss later in a simpler example. 

4-2. First order transitions 
4-2.1. Ordered models 

We now turn to the discussion of first order phase transitions. Perhaps the 
simplest mean field model without disorder that shows such a transition is a 
generalization of the quantum Curie- Weiss model to 3-spin interactions [209, 
103, 210]: 

T N N 

5 = -^E sf^aj-r^sf . (47) 

i.j,k=l i=l 

Like in the Curie- Weiss model the mean field nature of the model allows for 
an exact computation of its free energy density. A simple way to obtain the 
self-consistency equation on the order parameter (that can be formally justified 
via a path integral representation of the model [209, 103]) is to replace two a z 
by their average m and obtain a single site Hamiltonian H = — J m a" — Yo x , 
from which we get the mean field equation 

m = Jm tanh(^V^ 2 ™ 4 + T 2 ) . (48) 

v j 2 m 4 + r 2 

The paramagnetic solution m = of this equation always corresponds to a 
local minimum of the free energy. A ferromagnetic (m > 0) solution however 
appears discontinuously on a spinodal line in the (/?, T) phase diagram. The free 
energies of the two locally stable phases cross on a first order transition line, 
distinct from the spinodal of the ferromagnetic phase. The phase transition is 
here characterized by a jump of the magnetization. Correspondingly, the first 
derivative of the ground state energy with respect to Y has a jump, hence the 
transition is of first order. 
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In this case, one can show analytically that the gap is exponentially small at 
the phase transition, A min ~ cxp(— fJ,N), and compute analytically the coeffi- 
cient /i [209, 103]. The reason behind this is indeed quite simple. The transition 
is characterized by an (avoided) crossing between two different eigenvalues: the 
state \P) corresponding to the paramagnetic solution, which is the ground state 
at large V, and the state \F) corresponding to the ferromagnetic solution, which 
is the ground state at small T. One can see that these two states have an 
overlap (F\P) ~ exp(— fiN). Naturally, the matrix element of the Hamiltonian 
between these states is 7 = (F\H\P) ~ exp(— fiN), and an analysis of the two 
level system similar to the one of Sec. 2.3.3 leads to the exponential scaling of 
the gap at the phase transition. Therefore, a quantum annealing of this model 
requires an exponential time to find the ground state (despite the latter is again 
the trivial ferromagnetic one). It is also interesting to remark that performing 
a quantum annealing (as well as a classical annealing) of this model for a finite 
(with respect to N) time r always leads to a final energy that is extensively 
higher than the ground state one, even when r — > 00 (after N — > 00) [103]. This 
means that even the problem of approximating the ground state energy of this 
model through quantum annealing is very hard. The effect of spinodals in mean 
field models with first order phase transitions on quantum annealing has been 
recently discussed in [103]. 

Note that the physics of finite dimensional models undergoing a first order 
phase transition is quite different from the one of mean field models. This is 
because in finite dimension the dynamics around a first order transition is dom- 
inated by nucleation events that are absent in the mean field treatment [161]; 
this should affect the scaling of the gap, as well as boundary conditions in some 
unidimcnsional models [211]. We shall not discuss these issues further because 
in the following we will be mainly interested in mean field models. 

4-2.2. Disordered models 

We now turn to disordered spin glass models that show a quantum first 
order phase transition. The simplest such model is obtained by introducing 
disordered couplings in Eq. (47), which amounts to perform the same step that 
leads from the Curie- Weiss to the SK model, or to add a transverse field to the 
fully connected p-spin model of Eq. (20). The resulting Hamiltonian is the one 
of the 3-spin quantum spin glass [212, 35]: 



Here, the Jijk are quenched i.i.d. Gaussian variables, with zero average and 
Jfj k = J 2 /(2N 2 ); hence, the couplings can be positive or negative, leading 
to frustration. The classical (r = 0) thermodynamics of these models is very 
similar to the one of random optimization problems, as was discussed in Sec. 3.4, 
and is described by the so-called "Random First Order Transition" (RFOT) 
theory; for this reason they will be particularly relevant for the rest of the 
discussion. 




i<j<k 



(49) 
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Figure 9: (Left panel) Spectrum of the QREM as a function of T [38]. Red dots represent the 
results from exact diagonalization of a system with N = 20, dotted lines are the analytical 
values from lowest order perturbation theory. The inset shows the scaling of the minimal gap 
as a function of the size TV. (Right panel) Phase diagram as a function of temperature and 
transverse field [34, 38]. 



This model and similar ones were studied through techniques that combine 
replica and Suzuki- Trotter methods [212, 35, 37, 36], and it was shown that the 
model undergoes a first order quantum phase transition at low temperatures. 
In order to avoid introducing the replica method, that is not relevant for the 
present discussion, in the following, instead of discussing the Hamiltonian (49), 
we consider the simplest representative of the RFOT universality class, namely 
the Quantum Random Energy Model (QREM) [34, 38, 41]. This model is just 
the classical Random Energy Model (REM) introduced in Sec. 3.3 to which one 
adds a quantum transverse field, and it can be thought as a model similar to 
Eq. (49), but with interactions involving p spins in the limit p — > oo [140]. It is 
described by the Hamiltonian 

<t i 

where E(g_) are i.i.d. random variables, extracted from a Gaussian probability 
density with zero average and variance N/2. 

The complete phase diagram of the QREM as a function of T and T has 
been obtained in [34] by means of the replica method, and is reported in Fig. 9. 
The existence of a first order phase transition in (a slight variant of) this model 
has been rigorously confirmed in [41]. This phase diagram can be obtained via 
a series of extremely simple arguments [38]: 

• Extreme cases 

We start by discussing the two limiting cases (i) T = and (ii) T = oo. 

(i) When T = the model has a phase transition as a function of the 
temperature [140], that we already discussed in Sec. 3.3 and we briefly 
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recall here. The micro-canonical entropy density is s(e) = log 2 — e 2 on the 
interval of energy densities e = E/N where it is positive, i.e. [eo, — eo] with 
e = — Vl°g2. In consequence a condensation (or glass) transition occurs 
at the critical temperature T c = l/(2\/log2): at high temperatures the 
free energy density is /rem (T) = — -jy — T log 2 and an exponential number 
of configurations contribute to the partition function. On the contrary at 
low temperatures the Gibbs measure is concentrated on a finite number 
of configurations of energy density eo, and /rem {T) — eo = — ^/log2; the 
entropy density then vanishes. 

(ii) When T ^> 1 the system corresponds to TV independent spins in a 
transverse field, i.e. a simple Quantum Paramagnet (QP). As a conse- 
quence the free energy reads /qp(T, T) = — Tlog(2 cosh(F/T)) and the 
entropy density s(e) is the logarithm of a binomial distribution on e € 



• Perturbation theory 

We now discuss the perturbation theory around these two limiting cases. 
Consider first a classical configuration a with a negative energy density 
e(a) = E(a)/N < 0. The energy density e(a,T) of the corresponding 
cigenstate of H can be computed order by order in T, treating the trans- 
verse field as a perturbation, and this yields: 



as typical configurations that are reached by a single spin flip from a 
low energy configuration of the REM are the most numerous ones, with 
vanishing energy density (this will be further explained in Sec. 4.3.2). 

In the opposite limit r> 1, the eigcnstate of the pure transverse field arc 
degenerate so we must use degenerate perturbation theory. Consider the 
space of eigenstates where N — k spins are aligned in the transverse field 
direction. This space has degeneracy (^) . Here we consider the classical 
part of the Hamiltonian, Hp — ^ a E{a) \s)(s_\ as a perturbation. It is 
easy to show that the restriction of this Hamiltonian to each degenerate 
subspace is a random matrix, whose elements are Gaussian variables with 
zero mean and variance N/2 N . Therefore this is a random matrix belong- 
ing to the Gaussian Orthogonal Ensemble (GOE) and its spectrum is the 
usual semi-circle law with eigenvalues ~ N/2 N . We obtain that on the QP 
side, at first order, the energy levels are those of free spins in a transverse 
field, with exponentially small corrections. See [213] for the second order 
computation. 

The important outcome of these considerations is the vanishing as N — > oo 
of the perturbativc corrections around the two limits T = and r > 1. 
hence the energy, entropy and free energy densities are not modified. The 
free energy density of the QREM is thus /qrem = min[/ RE M(r), /qp(T, T)}, 



[-r,r]. 




(51) 
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which leads to a first order phase transition when the values of the free 
energies of the two competing phases become equal. The spectrum and 
the phase diagram of the model in the plane (T, T) are shown in Fig. 9. 

• Gap 

A good approximation of the minimum gap is given, as in the ordered 
case, by considering a two level problem similarly to (15) where the space 
is that spanned by the ground states of the classical part of H and of 
the transverse field, denoted respectively \Eq) (that corresponds to the 
classical state of minimal energ y) and \QP) = 2~ N / 2 £ CT \a). The diagonal 
matrix elements are the perturbed energies. The off-diagonal elements are 
proportional to the overlap (E$ \QP) = 2~ JV / 2 . From this we obtain that 
A m i„ oc 2~ N / 2 and the accuracy of this scaling is shown in the inset of 
Fig. 9. 

The phase diagram of the QREM shares strong analogies with the results 
obtained in other quantum spin glass models belonging to the RFOT class [212, 
35, 37, 36], consistently with the fact that the REM is a good approximation for 
more complex systems. In all these problems the classical glass transition occurs 
as a function of the temperature. At T = the system is in the glass phase 
and the classical ground state is not extensively degenerate. We stress here that 
this is a crucial difference with respect to many other optimization problems, 
like fc-SAT, where the glass transition also arises at T = as a function of the 
density of constraints. In these cases the entropy density is non-vanishing at 
zero temperature and the transition is an entropic phenomenon. We will show 
next how the role of entropy can be taken into account in a simple extension of 
the REM, the random subcubes model we already introduced in Sec. 3.6. 

Note the analogy of the computation of the gap in the QREM with the 
one of the Grover problem discussed in Sec. 2.3.5. However, in the QREM the 
classical intensive ground state energy has fluctuations of order l/y/~N, which 
induce similar fluctuations of the location in T of the minimal gap. Hence, in 
this case the optimal schedule discussed in Sec. 2.3.5 for the Grover problem 
(that is based on the exact knowledge of the location of the minimal gap) cannot 
be applied and a quadratic speedup is not achieved by the QAA in the QREM. 
Hence, the behavior of a QAA (as well as that of a classical annealing) is the 
same as in the ordered case: finding the ground state takes a time ~ 2 N , 
exactly as in an exhaustive search. An annealing over a finite time will lead to 
an extensive residual energy. 

4-3. Level crossings and localization on the hypercube 

In sections 4.1 and 4.2, we presented an overview of several simple models 
that show different phenomenologies: first and second order quantum phase 
transitions, associated to different scalings of the gap at the transition. More- 
over, we announced that some of these models are characterized by phases where 
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the gap is everywhere exponentially small due to an accumulation of level cross- 
ings [102, 39, 40]. In this section, we shall give a more detailed description of 
this level crossing phenomenon. 

The approach wc shall use here has its roots in the physics of Anderson lo- 
calization. In [139], Anderson considered a particle hopping on a d-dimcnsional 
cubic lattice of N = L d sites, and subject to a random disordered potential, as a 
good starting point for the comprehension of transport properties of metals and 
of the metal-insulator transition. The Anderson model describes an electron 
hopping in a disordered environment and the Hamiltonian reads 

N 

H AM = -* X! (cfo + cfii) + ^ ^% , (52) 
(i,f) i=1 

where the first sum runs over the edges of the lattice, i.e. pairs of sites at dis- 
tance 1, the second sum runs over all N sites, c\ is the creation operator at 
vertex i of the lattice, and £j arc i.i.d. random local energies, taken from a given 
distribution. Anderson showed that, depending on dimensionality and on the 
strength of the disorder, the eigenstatcs of the Hamiltonian (52) can be extended 
or localized in real space [139]. What is particularly interesting for the present 
discussion is that the spectral properties change completely in the extended or 
localized regions of the spectrum. Indeed, extended states are typically sepa- 
rated by much larger gaps than localized ones. In the literature on Anderson 
localization, this is sometimes referred to as level repulsion. Level repulsion is 
suppressed for localized states exactly because the matrix elements that con- 
nect them are much smaller, hence the 7 in Eq. (15) is much smaller leading 
to smaller gaps. Therefore, one might expect that the presence of avoided level 
crossings leading to exponentially small gaps could be interpreted as some kind 
of localization phenomenon [40] . We discuss this point of view in detail in the 
rest of this section. 

4-3.1. A different view on the QREM: the Anderson model on the hypercube 

A simple observation is that the transverse field operator Hq = — £\ of has 
non-zero matrix elements between two states \q) and \a') if and only if the Ising 
spins configurations a and er' differ on exactly one variable Oi . Considering the 
2 N configurations a E {— 1, +1} as the vertices of the N dimensional hyper- 
cube, and defining the Hamming distance d(a,a') as the number of different 
bits between the two configurations a and a' , one can rewrite any Hamiltonian 
of the form (50) as 

h = J2 e (<l)\<l)(<l\ - r X (1^) & + |£> <d\) . (53) 

where the second sum runs over pairs of neighboring configurations on the hy- 
percube, i.e. such that d(a,c/) = 1. In this formulation the QREM discussed 
in Sec. 4.2.2 is exactly the Anderson model on the hypercube, as the energies 
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E((j) are random i.i.d. variables, precisely as the of (52), the transverse field 
playing the role of the particle hopping term of the original Anderson model. 

As we discussed in Sec. 4.2.2, the energies E(a) of the QREM provide a 
disordered environment that induces localization on one of the vertices of the 
hypercube when the hopping T is not strong enough. This is shown by the 
fact that the low energy eigenstates at small enough V coincide with the classi- 
cal ones, hence with |er), at all orders in perturbation theory for N — ► oo, see 
Eq. (51) [38]. The crucial difference with the Anderson model is that here, the 
derealization transition coincides with the first order phase transition, and it 
happens via a level crossing between the localized and extended ground state. 
Moreover, in the localized (small L) phase, no level crossings are observed be- 
tween different states; this is clearly due to the fact that the energies of the 
lowest eigenstates do not depend on L, again due to Eq. (51) [38]. 

The analogy between the QREM and the Anderson model is strongly appeal- 
ing for the physicists community since it brings the field of quantum information 
and the one of Anderson localization of interacting systems in close contact [40] . 
On the other hand, the localization in the QREM is "extreme" in the sense that 
in the localized phase the eigenstates of the Hamiltonian coincide with the clas- 
sical states. This is due to the fact that the Hamiltonian is non-local and flipping 
a spin typically costs an extensive energy. 

Therefore, several interesting questions remain open. In local Hamiltonians 
(Sec. 2.2.2), one expect that for finite L, states are always delocalized over 
an exponential number of states of the computational basis (see the discussion 
in [214]). What is then the meaning of many-body localization? Is it possible to 
observe, in more general models, level crossings in the "localized" phase? Does 
derealization always happen through a first order transition? And finally, what 
happens when the classical ground state is exponentially degenerate, unlike in 
the QREM? In the following we try to answer some of these questions. 

4-3.2. A mechanism for level crossings between localized states 

Let us first consider the case of local Hamiltonians in the sense of Sec. 2.2.2. 
A mechanism that induces level crossings between localized states was proposed 
independently by Altshuler et al. [40] and by Amin and Choi [39]. It relies 
on the fact that in optimization problems with local interactions, the diagonal 
energies E(a) are not uncorrelated as in the QREM. Therefore, a careful choice 
of the classical energy function can lead to level crossings. We now review this 
construction. 

In [40] , a classical random energy function that is a sum of local interactions 
was considered (namely, the Exact Cover problem defined in Sec. 2.1.1). The 
analysis starts by choosing an instance of the problem with M — 1 clauses, such 
that there are at least two isolated (i.e. separated by an Hamming distance of 
order N) solutions ct x and a 2 of all the M — 1 clauses. These configurations 
represent degenerate eigenvectors for L = 0. However, as soon as V > the 
two ground state energies must split. Let us call l-E'i(r)) and \E- 2 (T)) the two 
eigenvectors that transform continuously into \q_i) and \g_ 2 ) for T — > 0, and 
Ei(T) and E2(T) the corresponding energies. 
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The splitting of the two solutions can be computed using perturbation theory 
for small enough T, as we already did for the QREM. One can write for any 
given non-degenerate classical eigenstate |cr) : 

oo 

E(T 1 a)=E(a) + ^ 2n Fn(2.) , (54) 

n=l 

with some coefficients F n (a). The effect of the presence of two degenerate so- 
lutions will appear only at order n oc TV if the two solutions g_ Y and g_ 2 have ex- 
tensive Hamming distance, and therefore it was neglected in the discussion [40]. 
Let's look to the first order as an example. It has the form 



1 

a 1 :d(cr,(T') — 1 



Fig. 10 highlights the crucial ingredient in the construction of Amin and Choi [39] 
and of Altshuler et al. [40], and compares it to the QREM. In the latter, the 
classical energies are uncorrelated and for low energy eigenstates, the coefficient 
F\ turns out to have a finite limit for TV — >• oo. This is because a low energy 
configuration is connected by a single spin flip to TV configurations (hence there 
arc TV terms in the sum), but those typically have an extensive energy differ- 
ence above it (hence the denominator is of order TV). Considering intensive 
energies, the correction is therefore of order T 2 /N as given in Eq. (51). The 
crucial difference for correlated energies E(a) that are sums of local terms is 
that a spin flip always leads to a finite energy difference with respect to the 
starting point. Therefore, the denominators in the perturbative expansion are 
of order 1 , and F\ ~ TV is of the same order of the classical energy. Extending 
the argument to higher orders one easily sees that all orders in perturbation 
theory are proportional to TV, if the energy is the sum of local interactions, and 
therefore contribute to the r-depcndence of the energy levels. 

Hence, recalling that both g_ x and g_ 2 arc assumed to have zero classical 
energy and calling F^ 2 = Fn^a^) — F n (a 2 ), we have: 

oo 

E 1 (r)-E 2 (r) = Y,r 2n Fn 2 ■ (56) 

n=l 

It was argued in [40] that for random problems the coefficients F n (a) have the 
same average over the disorder; hence, F^ 2 has zero mean and naturally one 
can assume that (F^ 2 ) 2 ~ TV. This leads, on average, to 

oo 

Bi(T) - E 2 (T) = VN^ 2n fn , (57) 

n=l 

where /* 2 are finite for large TV. Eq. (57) shows that, if the first non-zero 
coefficient is negative, one can find a small enough T* such that E 2 (]?*) — 
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Figure 10: A comparison between the classical energy function of the QREM (left panel) and 
the one of the correlated system studied by Amin and Choi [39] (right panel). The horizontal 
axis is a sketch for the 2^ dimensional space of classical configurations, while the vertical axis 
is the corresponding energy density. Each point represents a distinct classical configuration, 
and neighboring points have Hamming distance equal to 1. 



Now, we can add an M-th clause to the problem and fix AE in such a way 
that the new clause introduces at most a penalty AE on the classical energy (in 
the Exact Cover problem, AE = 4). Then, it still holds that E 2 {T*) > E^T*). 
At the same time, there is a finite probability that the additional clause will 
be satisfied by a 2 but not by so that Ei(Q) > £2(0). Because in the 
Hamiltonian there are no particular symmetries, the matrix element between 
|-Ei(r)) and \E2(T)) will be non-zero and the introduction of the additional 
clause induces an avoided level crossing, as in Eq. (15). 

Once again, the avoided level crossing is associated with an exponentially 
small gap because by assumption the two solutions have an Hamming distance 
of order N, hence they can only be connected at order N in perturbation theory, 
leading to a matrix element of order T N , i.e. exponentially small. It is important 
to stress that in this construction the crossing happens for T < T* ~ TV -1 /' 4 " \ 
hence at very small T in the thermodynamic limit. For this reason, these cross- 
ings have been called perturbative crossings in the literature [40, 41, 215]. This 
scaling was not clearly found in the numerical experiments, but this was at- 
tributed to the small exponent, visible only for large N [40]. 

A similar phenomenon was discussed by Amin and Choi [39]. Their con- 
struction is based on a system whose classical energy function has a deep but 
narrow energy minimum, and a secondary local minimum which is higher in 
energy but wider (see Fig. 10, right panel). This means that around the sec- 
ondary minimum, flipping a spin costs less energy. The denominators in the 
perturbation theory, Eq. (55), are smaller around the secondary minimum than 
around the global one. In turn, the secondary minimum lowers its energy more 
quickly under the action of the transverse field and eventually crosses the clas- 
sical minimum. 



£i(r*) > AE for any finite AE 



In fact, 



r* 



AE ^ 



(58) 
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The analysis of [40, 39] points out an important mechanism that can induce 
level crossings representing a serious bottleneck for quantum algorithms, that is 
a missing ingredient of the QREM, which is the simplest Anderson-like model 
of localization on the hypercubc: namely, the fact that in systems with local 
interactions the diagonal energies E(a) are correlated, and most importantly a 
single spin flip always lead to a finite change in the classical energy, which is 
not the case in the QREM. Thanks to this, the energy densities of the classical 
eigenstates have a non-trivial perturbativc expansion in T. Therefore, one can 
find particular realizations of the disorder, such that the energy of a classically 
excited state decreases faster, as a function of T, than the ground state energy, 
leading to a crossing at small T. 

This mechanism is for the moment only understood in perturbation theory 
(whose validity for these systems has been criticized [216]), and therefore it 
holds whenever perturbation theory holds, that is, if at small enough T the full 
eigenstates of the quantum problem remain close enough to the classical eigen- 
states. This is what has been called a "many-body localization" phenomenon 
in [40] . The other important ingredient is a very particular construction of the 
instances of the problem, that admit only two solutions. But, as we already dis- 
cussed in Sec. 3.9, typical instances of generic random optimization problems, 
even close to the satisfiability threshold, have an exponentially large number 
of solutions, and so we expect that non-degenerate perturbation theory should 
not hold, and the spectrum should be much more complex. Therefore we would 
like to understand what happens generically in problems that exhibit multiple 
and not necessarily isolated solutions. In particular, do the avoided crossings 
remain finite and isolated in T (hence leading to singularities in the ground 
state energy for N — > oo) or do they proliferate and accumulate, leading to a 
continuum of level crossings and a gaplcss phase? The latter question is par- 
ticularly important, because it has been argued that a finite number of level 
crossings can be eliminated by suitable redefinitions of the quantum Hamilto- 
nian [217, 218, 42, 219, 43, 215, 41]. We will address it in the next section. 

4-4- Level crossings and the role of entropy: the random subcubes model 

Motivated by the previous discussion, we now analyze a simple extension of 
the REM that takes into account the role of the massive ground state degen- 
eracy: the Random Subcubes Model (RSM) introduced in its classical version 
in [172]. Given that this model can be fully solved and reproduces most of the 
phenomenology we are interested in this section, we will discuss its properties 
in some detail. Some of these results have been published in [47]. 

The cost function (or problem Hamiltonian Hp) of the classical model has 
been defined in Sec. 3.6. We recall here that there are 2 Ar ( 1-a ) random clusters 
(that have the topology of sub-hypercubes of the total Hilbert space, hence the 
name subcubes); when these clusters are disjoint (the regime of interest here) 
configurations belonging to a cluster A have a random classical energy eo(A). 
Configurations that do not belong to any cluster have classical energy V and 
we assume that V ^> max^eo^). The classical properties of the model have 
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Figure 11: A comparison of the energy density landscape of the REM (left panel) and RSM 
(right panel), using the same conventions as in Fig. 10. The main difference is that in the 
REM, low energy configurations are isolated and typically surrounded by configurations that 
have extensively larger energy. Conversely, in the RSM low energy configurations arc arranged 
in clusters, each containing 2 Ns degenerate neighboring configurations. 

been discussed in Sec. 3.6, where it was shown that at large enough a > a scp 
the space of low energy configurations is decomposed into a set of disconnected 
clusters separated by large energy barriers. This is illustrated in Fig. 11 that 
should make the difference between the REM and RSM evident. All our analysis 
of the Quantum RSM (QRSM) will be restricted to the region a > a sep when 
the clusters are well disjoint (see Sec. 3.6). 

The main result of this section will be that quantum fluctuations, combined 
with the cluster structure, give rise to a scries of level crossings induced by a 
combined cnergctic-cntropic effect. Before going into the details, it is useful to 
give an overview and illustrate the way in which the model will be analyzed. In 
Sec. 4.4.1 we will discuss the spectrum of the clusters at finite TV in the limit V ~> 
oo. In this limit the Hamiltonian is block diagonal, each block corresponding to 
one cluster. The spectrum is characterized by true level crossings between states 
belonging to different clusters. The level crossings are due to the interplay of the 
classical energy and the classical entropy of the clusters. Quantum fluctuations, 
indeed, favor more entropic clusters. In Sec. 4.4.2 we will consider the case of 
finite V (still at finite N). A finite V reintroduces a lot of additional states 
that have to be taken into account. They allow to connect the clusters by single 
spin flips, therefore the Hamiltonian is no longer block diagonal. We will treat 
this situation by perturbation theory and variational arguments to show that 
a finite (large) V induces only minor modifications with respect to the infinite 
V case. In Sec. 4.4.3 we will investigate the low energy spectrum obtained by 
exact diagonalization for finite system sizes, and show that it is characterized 
by several avoided level crossings. In Sec. 4.4.4 and 4.4.5 we will consider the 
thermodynamic limit N — > oo. In Sec. 4.4.4 we will focus on the ground state, at 
T = 0. We will show that there is a first order phase transition that separates 
a Quantum Paramagnetic (QP) phase, at large T, from a Spin Glass (SG) 
phase, at smaller T, like in the QREM. In the spin glass phase the ground state 
continuously changes from one cluster to the other as a function of T, because of 
the level crossings between different clusters; the latter accumulate for TV — > oo 
giving rise to a unique SG phase, and are therefore distinct phenomena with 
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Figure 12: Schematic Hamiltonian matrix representing a finite size realization of the QRSM 
with 3 clusters (red, blue and purple components). The biggest green sector represents states 
that do not belong to S. The Hamiltonian has zero matrix elements between states belonging 
to different clusters because for a > ce Bep their Hamming distance is bigger than one. The size 
riA of each cluster block is fixed by its entropy = 2 Na ^ A ^ while the size of green component 
is much larger ny ~ 2 . We indicated with T the sectors of the Hamiltonian where there 
are non-zero off-diagonal matrix elements; still T connects only classical configurations at 
Hamming distance 1, therefore the matrix is very sparse in these blocks. 

respect to the first order transition. Finally in Sec. 4.4.5 we will consider the 
case T > 0: we will show in particular that quantum fluctuations promote the 
glass transition. In Sec. 4.4.6 we summarize and we comment the results. 

4-4-1- Spectrum of the cluster Hamiltonian 

Wc will now study the spectrum of the quantum Hamiltonian H = Hp+THg 
as a function of T, and from now on we focus on the region a > a scp where 
clusters are well separated (see Sec. 3.6), which is the most interesting for our 
purposes. The computation of the spectrum for a < a sep is more complicated, 
because in this region the clusters have overlaps and the arguments below do 
not apply straightforwardly (although they might be generalized for a > ay 
where the overlaps are exponentially small [172]). A schematic example of the 
Hamiltonian describing a finite system with three clusters in the regime where 
the clusters are well-separated is shown in Fig. 12. 

Remember that we call iS the set of all classical configurations that belong to 
at least one cluster (Ha being the Hamiltonian of a cluster), and have therefore 
classical energy extensively smaller than NV; configuration that do not belong 
to S have energy NV (and are described by the Hamiltonian Hy)- The total 
Hamiltonian is H = J2a Ha+H v +THq. We consider first the ("hard") V -> oo 
limit where Hp is infinite for the states that do not belong to S: then we can 
project out these states from the Hilbert space and look to the restriction of 
H = J^a Ha + FHq on 5, which contains 2 tot states. Because the matrix Hq 
only connects configurations at unit Hamming distance, and different clusters 
have distance of order TV, the Hamiltonian H has no matrix elements connecting 
different clusters. Therefore we can diagonalize H separately in each cluster. 
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Figure 13: Low energy spectrum for a system with N = 15 and 3 clusters at Hamming distance 
larger than 1 such that {{s(A t ), e(Ai)) i=1 , 2 ,3} = {(2/15, 0); (4/15, 0.1); (6/15, 0.3)}. 
(Left panel) Spectrum in the V = oo case. To each cluster Ai corresponds the spectrum of 
Ns(Ai) free spins in a magnetic field. 

(Flight panel) Partial spectrum for finite V = 2 obtained by exact diagonalization, with a zoom 
on the ground state in the inset. Green lines are the results of second order perturbation theory, 
the blue line is eQp = V — V. At small V the low energy spectrum is in good agreement with 
that at V = oo. For larger values of T, avoided level crossings appear (marked by red arrows). 
The crossings at smallest T involve different clusters. The largest crossing, instead, involves 
the ground state of the spectrum connected to the classical low energy spectrum (set 5) and 
the ground state of the V-band; this crossing becomes a true first order phase transition in 
the N — > oo limit. 

The restriction of H to a given cluster A with Ns(A) free spins is equal to 
Ha plus the Hamiltonian of Ns(A) uncoupled spins in a transverse field, its 
spectrum is hence made of levels 

E k (A) = Ne (A) + (2k - Ns(A))T , k = 0, • • ■ , Ns(A) , (59) 

each ( Ns f. A ^) times degenerate. In particular the lowest level has energy per spin 
eas{A) = eo(A) — Ts(A), therefore the energy of clusters with larger entropy 
decreases faster with T. In this regime then one expects level crossings between 
states belonging to different clusters. In the situation where bigger clusters at 
r = have larger classical energy which is the case for most random optimiza- 
tion problems, the level crossings concern the ground state and at T = each 
crossing corresponds to a global rearrangement of the system. A simple example 
of a spectrum in the V = oo limit regime for a finite system containing three 
clusters is shown in the left panel of Fig. 13. Note that as long as the clusters 
are well separated, due to the V — > oo limit, there are no corrections in the size 
of the system. The crossings are not avoided and the degeneracy of the states 
is not removed, due to the complete independence of the Hamiltonian sectors 
describing each cluster. 

4-4-%- Quantum paramagnetic state 

Next, we consider a "soft" version of the model in which V is finite (still 
with V ^> max j 4eo(A)). Therefore now H is defined on the full Hilbcrt space 
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%. Iii this case, in addition to the 2 Nstot energy levels discussed above (that we 
shall refer to as the 5-band), there exists another set of 2 N — 2 Nstot ~ 2 N levels 
(the F-band), whose energy is expected to be of order V at small V. 

For the states in the 5-band we use perturbation theory in T. As soon as the 
transverse field is switched on a first order correction in T to the states in the S- 
band is present. This correction comes from the partial lifting of the degeneracy 
within the cluster and it is given by the spectrum Ek(A) in Eq. (59). A second 
order correction is induced by the presence of the V-band at finite V. In order 
to compute it one can apply perturbation theory assuming as unperturbed basis 
the one that diagonalizes the perturbation Hq inside each cluster. In particular 
we are interested in the correction to the lowest energy level ecs{A) in the 
clusters whose state \GS{A)) is given by all free spin polarized in the direction 
of the field. Then the correction is 

aeT i = v MH Q \GS(A))\i = T 2 (l -s(A))N 

A \ma E ^- E os(A) ~NV-Ne GS (A) 

and at any finite order n the correction to the energy per spin is 0((T 2 / (NV)) n ), 
so it vanishes in the thermodynamic limit. This mechanism is similar to the 
QREM and is due to the fact that states at the boundary of the clusters have 
extensive larger energy. 

To study the lowest energy level in the F-band e GS( y) it is convenient to 
rewrite the Hamiltonian in the following way: 

H = NVT-rY,o?-N , £,(y-e (A))\A){A\=H Q p-H s , (61) 



H QP H s 

where / is the identity and |^4)(A| = X^eA indicates the projector over 

the cluster A. Hqp acts on the entire Hilbert space while Hs acts only on 
the subspace spanned by the clusters. This form aims to interpret Hs as a 
"perturbation" over Hqp which describes a system of N free spins in a trans- 
verse field (with a shift JVV" in the energy). However the "perturbation" is 
not in the strength of the energy, which may be large, but in the number of 
states that are involved. Note, in fact, that Rank(iTs)<CRank(ifgp), being 
Rank(^ 5 )= K = 2 Na * ot and Rank(H QP )= 2 N . This, together with the fact 
that the perturbation matrix is positive defined (it shifts some states all in the 
same direction) allows to apply the results of small rank perturbation analy- 
sis [220] in order to study eos(v)- From these results we can safely say that 

E k Q - p n < E% < E k QP forfc = l,...,2 iV , (62) 

where Eq P and E H are respectively the fc-th eigenvalues of Hqp and H, and 
we assume Eq P = — oo when k < 0. In particular when T is small, eQs^ is 
larger than all the energies in the 5-band. This implies that 

V - r < e GS{V ) ■ (63) 
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The results from small rank perturbation (62) also shows that the spectrum of 
the F-band is close to the one of N free spins in transverse field with classical 
energy NV: 

Efy = NV + (2k - N)T , k = 0, • • • , N , 
with degeneracy close but not equal to ( fe ). We expect that the unperturbed 

ground state of Hqp, \QP) = 2~ N / 2 ^2 a \a) describes well the lowest energy 
level of the V-band &GS(v) and remains unaffected by the presence of the states 
in S for all T except from the region where it crosses the spectrum of S. The rea- 
son for this comes from the intuition that in absence of H$ the spectrum of Hqp 
is highly degenerate, especially in the middle of the band. Then, also comforted 
by the results of exact diagonalization, we expect that the states that recombinc 
the most in order to create the S-band when Hs is applied, are those belonging 
to the more degenerate part of the spectrum. On the contrary, \QP) is made of 
all spins aligned along T without degeneracy and thus it is weakly perturbed by 
H s . A rigorous study of this energy level is not possible, but we can use a vari- 
ational argument to understand its behavior. The state \QP) has exponentially 
small overlap with any state in the 5-band (ip(A)\QP) ~ 0(2- Na WI 2 ) and thus 
it gives an expectation value of if equal to (QP\H\QP) = N(V — T) + 0(2~~< N ) 
for some 7. If we interpret this as a variational upper bound on the true ground 
state of the F-band we get: 

e GS (v) <V-r. (64) 

Combining (63) and (64) we obtain 

e GS(v) =V-r + 0(2-^ N ) , (65) 

and the corresponding eigenvector remains up to exponentially small corrections 
the same \QP), which is uniformly extended in the basis \a). 

4-4-3- Exact diagonalization results 

We checked these predictions for the spectrum by means of exact diagonal- 
izations for a system made of N = 15 spins. The results are shown in Fig. 13, in 
the right panel. There we have plotted the spectrum of a system made by three 
clusters characterized by classical energy and entropy {(s(Ai), e(Ai))i = i t 2,3} = 
{(2/15,0); (4/15,0.1); (6/15,0.3)} and V = 2. The plot shows that for small T 
the states in the T^-band do not affect those in the set S, whose spectrum is in 
good agreement with that at V = 00, in the left panel. At larger T, avoided level 
crossings first appear between the ground states of different clusters. Finally, an 
avoided crossing happens with the ground state of the l/-band, whose slope in T 
is much larger due to the big entropy which characterizes this sector. We have 
also plotted in green the analytical result that we obtain up to second order in 
perturbation theory for the lowest energy level of each cluster and in blue the 
energy of the quantum paramagnetic state. We see that the true ground state, 
crossing after crossing, well interpolates between all these curves. Since the clus- 
ters have Hamming distance proportional to N, we expect all these crossings to 
be avoided at finite N producing exponentially small gaps [39, 40, 41]. 
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Figure 14: Results for the QRSM in the region ct Bcp < a < a B . As an example we choose 
(following [172]) p = 0.7, a = 0.85, and g(e ) = [2 + e /e m - (e /e m ) log(e /e m )]/3 for 
eo S [0, e m ] with e m = 0.1. 

(Le/i panel) Energy of the SG ground state [Eq. (66), full line] and of the QP state eQp = 
V — r for V = 1 (dashed line). A first order transition between the two states happens at 
r ~ 2. (Inset) Level crossings in the SG state. For better readability we plot ego +rs max (0) 
[Eq. (66), full line] and show the energy eo — r[s ma x(eo) — Smax(0)] of two different clusters 
with eo = 0.05,0.2 (dot-dashed lines). 

(Right panel) Transverse magnetization function of T for the same parameters. The 

first order phase transition between the SG and the QP is manifested by a jump in m x , shown 
by a vertical dotted line. The value of m x for the SG is the solid black line, while that of the 
QP is the dashed blue line. Note that the latter is bigger because it corresponds to a more 
cntropic phase. 



4-4-4- Level crossings in the thermodynamic limit 

We discuss now the zero temperature phase diagram of the model for a > 
a scp and N — > oo. Following Sec. 3.6, to get a meaningful thermodynamic 
limit, the number of clusters of energy eo is set to 2 7V ( 1 ~ Q ) £, ( e °), where gifio) is 
an arbitrary increasing function of eo G [0, e m ] (as in most random optimization 
problems). We assume that g(e m ) = 1 so the total number of clusters in S 
is still 2 Nt - 1 ~ a \ As discussed in Sec. 3.6, the complexity of clusters of energy 
eo and entropy s is S(eo,s) = (1 — a)g(eo) — D(s\\l — p), and it vanishes at 
Smax(eo) which is also an increasing function of eo. The 5-band, or spin glass 
(SG), ground state energy is 



The minimum is in eo = as long as T < T\ c = 1/(8^(0)). Above this 
value, the minimum is in a different eo for each value of T: in this region the 
ground state changes abruptly from one cluster to another upon changing T by 
an infinitesimal amount, similarly to what is called temperature chaos in spin 
glasses [184, 186]. Note that in some relevant cases the slope of g(eo) in eo = 
is infinite, therefore Ti c = and level crossings happen at all V. 

The energy eQp crosses the SG ground state given by Eq. (66), giving rise to a 



e S c = mm 

e S[0,e m ] 




(eo - Ts) 



(66) 



mm 

e e[0,e m ] 
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a 

Figure 15: Phase diagram of the model for p = 0.7, g(eo) as in Fig. 14, and P = 1 (full lines). 
The vertical line corresponds to a scp = 0.797 for this value of p. The higher T line is the 
first order transition between SG and QP. Above the lower T line a c (T, P = 1) the system is 
in the condensed phase. The condensation transition lines a c (T,P) are also reported (dashed 
lines) for different values of /3, showing that the non-condensed phase disappears for /3 — > oo. 
The complexity of the zero-energy clusters is (1 — a)g(eo = 0) = 2(1 — a)/3, hence one has 
a c (r = 0, P = co) = 4^ + | log 2 (2 -p) = 0.875. 



first order phase transition between the SG and the QP [38, 34, 35, 36, 37, 48] at 
a critical r cx V. As a consequence, the transverse magnetization m x = de/dT 
has a jump at the transition [48] (see the right panel of Fig. 14). Note that 
m x = s, thus the transverse magnetization is determined by the entropy of the 
ground state, and the entropy of the ^/-component is much larger than those of 
the clusters. 

4-4-5- Finite temperature: the condensation transition 

The previous analysis shows that in the region a sep < a < a c the pertur- 
bation THq has a dramatic effect. At T = 0, most of the states in S belong 
to one of exponentially many small clusters, while at any T > the few largest 
clusters of entropy s max have the smallest energy This is related to the fact 
that the presence of a transverse field introduces a correction to the energy 
that favors the more entropic clusters. A more complete picture is obtained 
by studying the model at finite temperature (recall that the classical model at 
finite temperature was studied in Sec. 3.6.3). It is convenient to separate the 
contribution of the two parts of the spectrum (the 5-band corresponding to 
the SG phase and the F-band corresponding to the QP phase) to the partition 
function, Z = Tr e^ 5 = Z SG + Z QP , with c = 2cosh(/?r): 

Zqp-^V^ =e-P NV c N , 

k 

Z SG ~Y. e ~ PEhiA) = I de ds2 N ^ e °> s h-P Ne °c Ns . 

A,k J 

Of course, Z^g reduces to the classical partition function in Eq. (39) for T = 0. 
The free energy is /qrsm = - (T/7V) log Z = min{/sG, /qp}, analogously to 
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what was found in [38] for the QREM (see the discussion in Sec. 4.2.2), with 
/qp = V - T log c and 

f SG = -T max [S(e 0j s) log 2 - /3e + s log c] . (67) 

e 6[0,e m ] 
s6[smin(eo),s max (eo)] 

The first order transition happens when the free energies /sg and /qp cross, 
while the condensation transition a c (T, T) happens when the maximum in Eq. (67) 
is attained in s max for the first time. In Fig. 15 we plot the lines a c (T, T) versus 
T for several temperatures. We observe that in the limit /? — > oo, the lines 
a c (T, r) shrink to the horizontal axis and the system is in the condensed phase 
for any T > 0. The first order transition to the QP phase happens for larger 
values of T at fixed temperature, and it is reported in the plot for (3 = 1. 

4- 4- 6- Summary 

Before presenting a more general perspective on random optimization prob- 
lems, let us summarize the results of this section. We introduced a simple 
toy model of a quantum optimization problem, the QRSM based on the RSM 
of [172]. In the classical case r = 0, the model captures the essential structure 
of the space of solution of random optimization problems, and displays several 
phase transitions that are present also in more realistic problems such as fc-SAT, 
at least at large k. We explored the consequences of this complex structure on 
the spectrum of the quantum Hamiltonian at T > 0, and we showed that: (i) 
Quantum fluctuations lower the energy of a cluster proportionally to its size, 
(ii) Because the energy and the entropy vary from cluster to cluster, level cross- 
ing between different clusters are induced as a function of T in the SG phase, 
due to a competition between energetic and entropic effects. These crossings 
accumulate for TV — > oo in a continuous range of T, giving rise to a complex SG 
phase characterized by a continuously changing ground state and an everywhere 
exponentially small gap. (iii) At large r ~ V the SG phase undergoes a first or- 
der transition towards a QP phase, corresponding to the complete derealization 
of the ground state in the computational basis \q_). (iv) At finite temperature, 
there is a line of condensation transitions a c (r) that shrinks to T = at low 
temperatures: indeed, at zero temperature the condensation transition becomes 
abrupt. While at T = the space of solutions is dominated by an exponen- 
tial number of clusters of intermediate size, for any T > the biggest clusters 
contain the ground states. 

4-5. Phase transitions in quantum optimization problems: an attempt towards 
a general perspective 

Overall, the discussion of the previous sections shows that the low energy 
spectrum of quantum optimization problems can be very complex, and char- 
acterized by different level crossings: internal level crossings in the SG phase, 
or the crossing between the SG and the QP giving rise to a first order phase 
transition. Moreover, both entropic and energetic effects arc important. 
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Figure 16: Pictorial energy landscape of the REM (upper left panel), RSM (upper right 
panel), XORSAT on a satisfiable random regular graph (lower left panel) and fc-SAT (lower 
right panel), using the same conventions as in Fig. 10. 

Yet the previous discussion was based on a series of toy models (such as 
the QREM or the QRSM) or on the analysis of extremely simplified instances 
of random optimization problems. These problems were basically constructed 
ad hoc to exhibit the desired phenomenology The next task is therefore to 
demonstrate that these phenomena indeed happen in typical instances of realistic 
random optimization problems, such as those defined in Sec. 2.1.1. This will be 
the subject of Sec. 6, but it requires the introduction of sophisticated quantum 
statistical mechanics tools that we will discuss in Sec. 5. Before proceeding, 
in this section we want to complete the picture by presenting coherently what 
are the expected properties of the spectrum of generic random optimization 
problems. 

As was explained in Sec. 3, the classical energy landscapes of several random 
optimization problems have been recently characterized in much detail, thanks 
to important developments in the analysis of classical spin glasses [32, 129, 152, 
125, 136, 137]. These studies show that the classical energy is characterized by 
many "valleys" at the bottom of which local minima are found; we are now able 
to obtain a quite detailed quantitative characterization of the shape of these 
valleys [188]. We can use as running examples the random regular XORSAT 
problem, in the satisfiable phase where solutions exist, which is a representative 
of the class of locked models discussed in Sec. 3.9, and the random fc-SAT 
problem, which instead displays (for k > 3) all the transitions discussed in 
Sec. 3.7. In Fig. 16 we sketch pictorially the energy per spin as a function of the 
configuration for those models, and compare with the previously investigated 
toy models, the REM and RSM. The picture shows some very general aspects 
of mean field spin glasses, that are shared by all the models considered here: 

a. The energy landscape contains many local minima. 

b. The distance between low energy configurations belonging to the basin of 
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attraction of different minima is O(N). 

c. The height of the energy barriers separating two different minima is O(N). 

At the same time, there are some crucial features that are model dependent also 
within mean field models: 

d. The number of configurations around a given local minimum might be 
exponentially large (the entropy is positive) or not (the entropy is zero). 

e. The (intensive) energy change associated to a spin flip starting from a low 
energy configurations can be either 0(1) or 0(1/N). The latter case is 
the rule for Hamiltonians that are the sum of local terms, and in this case 
the "steepness" of the energy around a local minimum can depend on the 
minimum itself. 

Based on the previous analysis of toy models, we expect that each low energy 
cluster of the classical energy function E(g_) gives rise, under the action of a 
quantum term like a transverse field, to a set of states (whose size roughly 
corresponds to the classical entropy of the cluster), with an energy density of 
the form 

e(r) = e+ rO(l) + r 2 C(f) .... (68) 

Entropic effects Energetic effects 

The coefficients of both terms depend on the shape of the classical energy around 
the local minimum. It was shown explicitly for the QRSM that the coefficient 
of the linear term is directly proportional to the intensive entropy of the cluster. 
Therefore, if local minima are exponentially degenerate, this coefficient is 0(1) 
with respect to N, and its fluctuations from cluster to cluster are also 0(1). 
Moreover, it was shown in Sec. 4.3.2 that the coefficient of the quadratic term 
depends on the neighborhood of the local minimum. Fluctuations of the latter 
coefficient among different minima are 0(1/V r /V) in the example of [40] but 
might be 0(1) in other models. 

The cluster to cluster fluctuations of the quantum corrections will generically 
lead to avoided level crossings. Because of the huge number of different clus- 
ters, it is reasonable to expect that these crossings will accumulate for N — > oo 
leading to an everywhere gapless spin glass phase, as it was shown explicitly 
for the QRSM. This phase will cease to exist at large enough T, when the com- 
pletely delocalized state, corresponding to the quantum paramagnetic phase, 
will cross the cluster ground state. Generically we expect this crossing to be- 
come a first order transition in the thermodynamic limit (as in the QREM and 
in the QRSM), but the transition might also be of higher order depending on 
the model. We will discuss concrete examples in Sec. 6. In summary, we expect 
generic problems to display a complex spin glass phase for small T, separated 
by a quantum critical point from a simple quantum paramagnetic phase. 
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5. Methods 



The aim of this section is to review various methods that can be used to 
investigate quantum spin glass models. We will give a particular emphasis to 
the methods that are most efficient to treat the quantum versions of the family 
of random optimization problems introduced in Sec. 2.1.1 and Sec. 3.1. We 
will use these methods in Sec. 6 to obtain results on specific models such as 
the XORSAT and coloring problems. This section is organized as follows: we 
will start in Sec. 5.1 by introducing the classical cavity method, a framework 
that has been developed to study random ensembles of optimization problems 
and which has led to the understanding of their complex phenomenology, that 
we explained in Sec. 3.7. The next three sections (5.2, 5.3 and 5.4) will be de- 
voted to different approaches to the generalization of the classical cavity method 
to quantum models. Finally in sections 5.5 and 5.6 we will give some details 
on some more standard numerical methods, such as exact diagonalization and 
quantum Monte Carlo, that have also been used to obtain important informa- 
tions on these problems. Before entering in the core of the discussion let us give 
a more detailed overview of the rest of this section. 

As already explained in Sec. 3.2, disordered mean field models [32] can be 
roughly classified in two categories: fully connected ones, where each degree 
of freedom interacts weakly with all others, and finitely connected ones, with 
a finite number of strong interactions for each degree of freedom. The replica 
method has been originally devised for the former family, most notably for the 
Sherrington-Kirkpatrick model [142]. Its extension to the finite connectivity 
case [151, 129] has been more conveniently reformulated in terms of the cav- 
ity method [152], and applied in particular to random ensemble of Constraint 
Satisfaction Problems (CSP) [126, 130]. The classical cavity method is by now 
a well established technique, with many presentations in original research pa- 
pers [152, 221] and in a textbook [137], and rigorous proofs of validity in some 
cases [131, 132, 222]. For the sake of completeness, in Sec. 5.1 we provide a 
quick survey of the classical cavity method, before turning to the specificities 
of its quantum version. In this section we also discuss the general definition of 
random graph models and the key concept of replica symmetry breaking. 

Several quantum extensions of the cavity method have been recently pro- 
posed in a series of papers, able to treat for instance spin 1/2 models in pres- 
ence of a transverse field. Roughly speaking, these methods can be divided 
in three groups. Path Integral Quantum Cavity (PIQC) methods exploit a 
path integral representation in order to map the quantum problem into a clas- 
sical one and then make use of the classical cavity method [223, 224, 225, 48]. 
Operator Quantum Cavity (OQC) methods work directly with quantum op- 
erators [226, 227, 228, 229, 230]. Finally, Variational Quantum Cavity (VQC) 
methods propose a variational ansatz for the ground state wavefunction that can 
be represented in terms of a set of local parameters, and then use the classical 
cavity method to optimize the energy of the variational state [231]. 

PIQC is at the moment the only analytical method that was used to obtain 
concrete results on one of the random CSP defined in Sec. 2.1.1, namely XOR- 
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SAT in presence of a transverse field [48]. The goal of Sec. 5.2 is to explain the 
technical details of the PIQC at the level of one step of replica symmetry break- 
ing, that is needed for the solution of these problems [48], and to generalize it 
to arbitrary discrete quantum degrees of freedom 3 . We will rely on this method 
in Sec. 6 to present detailed results on the XORSAT model [48] and original 
results on the coloring problem. 

The main drawbacks of PIQC are that i) as any path integral sampling 
method, it is restricted to Hamiltonians that are not plagued by the "sign prob- 
lem" (or in other words that admit a path integral representation where the tra- 
jectories have positive weights), ii) it does not allow to work exactly at T = 0, 
but only to perform an extrapolation to T — > from finite temperatures and in) 
the resulting functional equations have to be solved using a statistical represen- 
tation, that is affected by fluctuations and/or finite size (of the representation) 
effects. These drawbacks could in principle be overcome by working directly 
with operators. In Sec. 5.3 we will describe several attempts to construct OQC 
methods [226, 227, 228, 229, 230]. In Sec. 5.4 wc will describe VQC methods 
that are specifically designed to obtain direct information at T — [231]. Al- 
though these attempts are extremely promising and already allowed to obtain 
very interesting results for simple ordered and disordered models, it turns out 
that for the moment they are too computationally demanding to be used to 
solve the problems of interest here (e.g. the XORSAT problem). This will be 
discussed in Sec. 6. 

5.1. The classical cavity method 
5.1.1. Factor graph models 

Let us recall some definitions of Sec. 2.1.1 and put them in a more general 
context. We consider a model with N degrees of freedom er, taking values in 
a finite alphabet X, for instance X = {— 1,+1} for Ising spins. We denote 
a = (o"i, . . . , erjv) the global configuration, and for a subset S C {1, . . . , N} we 
write a s = {o~i\i G S} the configuration of those variables. The model is further 
defined by its energy function E(a), that contains M interactions labeled by 
a = l,...,M: 

M 

E(a) = J2e a (a Sa ) , (69) 

a=l 

For each of the interactions da denotes the set of variables that interact through 
a with the energy term e a ; a convenient representation of such an energy function 
is provided by factor graphs [233], see left panel of Fig. 17. Each variable i is 
associated to a circle vertex (variable node), while interactions a are represented 
by squares (function nodes), an edge being present between interaction a and 
variable i if and only if the value of e a depends on Oj, i.e. if and only if i € da. 



3 Note that the PIQC has also been used in a condensed matter context, namely to inves- 
tigate quantum glassy phases of disordered interacting bosons on random lattices; a detailed 
explanation of the method in this case can be found in [232, 198]. 
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Figure 17: (Left panel) An example of a factor graph. (Right panel) Graphical representation 
of Eq. (71). 

We shall also use the notation di for the set of interaction nodes linked to the 
variable i, and call graph distance between two variables i and j the minimal 
number of interactions in a path along the factor graph between i and j. 

The Gibbs-Boltzmann measure at inverse temperature (3 of the model can 
be written as 

M N 

z = e n^fe^n^) > ( ? °) 

a£X N a=l i=l 

where the partition function Z ensures the normalization of the probability 
law, and where the interaction weights are given by w a = e~^ Ea . For future 
convenience we shall treat a slightly generalized case with weights Vi on the 
variables. For Ising spins, the latter can be thought as originating from local 
magnetic fields. 

Let us assume momentarily that the factor graph representing the model 
under study is a tree. Then the problem of characterizing the measure (70) and 
computing the associated partition function Z can be solved exactly in a simple, 
recursive way. One introduces on each directed edge i — > a from one variable 
i to an adjacent function node a a "message" jUi-> a , which is a probability 
measure on the alphabet X, that would be the marginal probability of <7j if 
the interaction a were removed from the graph. These messages are easily seen 
to obey the recursive (so-called Belief Propagation) equations depicted on the 
right panel of Fig. 17, 

Mi^a(o-i) = Vi(<Ti) Y[ ^w b (a db ) Y[ H-*b{vj) , (71) 

bedi\a !Lg b \i j€Ldb\i 

with Z 

i—ta ensuring the normalization of the law /i^— >a ■ 
On a tree factor graph there exists a single solution of these equations, which 
is easily determined starting from the leaves of the graph (for which the empty 
product above is conventionally equal to 1) and sweeping towards the inside of 



^ M N 
a— 1 i—1 
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the graph. Once the messages Hi^. a have been determined all local averages 
with respect to \x can be computed, for instance the marginal law for a single 
variable and for the set of variables around one interaction read respectively 



j€da\i 



with and z a defined by normalization. Moreover the partition function can 
be computed as 



N 



A I 



logZ = ^bg«i-5^(|aa|-l)logz a 

i=l 
iV 

= E l0 

i=l 
A/ 

- £(|0 a |-l)log 



aedig_ 



da\i 



(73) 



a=l 



5.1.2. Random ensembles 

In the definition of the energy function (69) and of the associated Gibbs- 
Boltzmann measure (70) we considered a single realization of the model. We 
turn now to random ensembles of such models. Their definition involves a prob- 
ability law on the integers, that we shall denote qa', for simplicity of notation 
we assume that all interactions a involve the same number k of variables. The 
factor graphs are then supposed to be chosen uniformly at random among those 
with TV variables and such that a fraction of variables is involved in d interac- 
tions. The number of interactions is then M = aN, with a related to the degree 
distribution q^ by the relation ak — ^2 d dq q . This definition is an hypcrgraph 
generalization of the random graph models with prescribed degree distribution, 
studied for instance in [234]. The Erdos-Renyi construction described in Sec. 3, 
where the M interactions are chosen uniformly at random among the (^f) possi- 
ble ones, is essentially equivalent to the one just described, with qd the Poisson 
distribution of mean ak. In the definition of the random energy function we 
also assume that the interaction and variable weights w a and Vi are chosen 
independently at random from two probability laws. We will denote E[-] the 
average with respect to the whole construction; the main objective in the field 
of disordered systems is the computation of the average free energy density / 
in the thermodynamic (large size) N — > oo limit, 



-Pf= lim ^E[logZ] 

N— yao iv 



(74) 
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as the free energy is a self-averaging quantity, which yields all relevant thermo- 
dynamic quantities by suitable derivatives with respect to its parameters. 

The cavity method allows to compute the thermodynamic limit of the free 
energy density for models defined on random graphs. It exploits some prop- 
erties of these random graphs, in particular their tree-like character: in the 
thermodynamic limit these random graphs converge locally to random trees. 
This means that the neighborhood within a fixed graph distance t of an ar- 
bitrary variable node i is, with a probability going to one as N — > oo with t 
fixed, a tree. The latter can be described by qd and an associated distribution 
Qd = (d + l)<7d+i/afc. Indeed the reference node i will have d interaction nodes 
around it with probability qd- Then each of the (fc — l)d variables at distance 
1 from i will have a number of descendants drawn independently from the law 
qd, and so on and so forth until t generations of vertices have been drawn (all 
degrees are drawn from qd except the first one with qd). One can notice that qd 
is the degree distribution of a variable node chosen uniformly at random, while 
qd corresponds to the selection procedure where one edge of the factor graph is 
chosen at random, say between i and a; then i belongs to d interaction besides 
a with probability qd- 

We explained above how statistical mechanics models defined on trees could 
be solved recursively. On the other hand we have just recalled that random 
graphs are locally tree-like; the cavity method is a set of prescriptions on how to 
exploit the local properties of the random graphs to make global predictions on 
the free energy density. Depending on the models, in particular on the amount 
of frustration between the variables they induce, different level of sophistications 
of the cavity method are necessary. 

5.1.3. Replica symmetric cavity method 

The simpler one, exact for instance for ferromagnetic, unfrustrated, models, 
goes under the name of Replica Symmetric (RS). In that case there exists a 
single pure state (or a small number of them simply related by explicit symme- 
tries of the model) in the Gibbs measure, which enjoys in consequence spatial 
correlation decay properties. The effect of the long loops which are present 
in the random graphs is then negligible in the thermodynamic limit, only pro- 
viding a self-consistent boundary condition. The recursive equations (71) valid 
for a single tree factor graph are given a probabilistic meaning to describe the 
thermodynamic limit of random graphs. More precisely, the order parameter 
of the RS cavity method is V(rf), the probability (over the disorder) that the 
messages /i,^ n in Eq. (71) (which are themselves probability distributions on 
X) are equal to rj. V obeys a self-consistent functional equation, which is more 
simply written in a distributional form as 

v = g(m,i, ■ ■ ■ ,m,k-i, ■ ■ ■ ,Vd,i, ■ ■ ■ ,vd,k-i,v,wi, . . . ,w d ) . (75) 

In this equation all the r)'s are drawn independently from V , and = denotes 
the equality in distribution between random variables. Moreover d is drawn 



78 



according to the law qd, the v and w a 's are independent copies of the variable 
and interaction random weights, and the function g in the r.h.s. is defined by 



r/(cr) 



1 



z({Va,i},V,{Wa}) 



v(cr) X 



[]la, t K,) J w a (a, <T ,1, ■ • • ,CTa,fc-l) , (76) 



r CT .ii6[l,*-l] \ o,t 



a=l 



^({^a.i}, {w a }) ensuring the normalization of rj. The RS prediction for the 
free energy in the thermodynamic limit is then 

~Pf = J"^ ^ E [l°g Z] = E lo g (z v ({Va,i} a % k d J 1] ,v,{w a } ae[hd] f) 

- a(k-l)E[log(zt({ Vi } ie[hk] ,w))] , (77) 

where the expectation is with respect the independent choices of d (with the 
degree distribution qa), of the 77's (according to the law V) and of the random 
weights v and iy a 's, and where the contribution of the variable and function 
nodes are defined by 



Zv({VaA' V A w a}) = 



I i«£[l, k-1] 

~[r]a,i( cr o;i) II w a{v,Pa,l, ■ ■ ■ ,cr a ,k-l) , (78) 
a , i la 

l~[[n l (crt)\w(<T 1 ,... 1 a k ). 

,—,(Tk \ i / 



One can also obtain the disorder average of local observables, for instance 
the average of the marginal probability for a single variable 0{ reads 



v { a ) HVa,i( a a,i) Y\w a {c,Cra,l, ■ ■ ■ ,<T a ,k-l) 



,ie[i,fc-i] 



J2 V ( a ') HVaA^aA n w a(0-',cr a ,i, . . . ,<J a .k-l) 



Other thermodynamic observables can be obtained in a similar fashion; one 
way to derive their expressions is to observe that the expression (77) for the free 
energy is variational, in the sense that the stationary conditions with respect 
to V coincide with the self-consistent condition of Eq. (75). In consequence one 
can use the partition function as a generating function of the observables to be 
computed, and take only explicit derivatives with respect to their conjugated 
fields in (77). 
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5.1. 4- Replica symmetry breaking 

The assumption of correlation decay that underlies the RS cavity method 
can fail in presence of frustration, for instance in the case of random CSP with 
a > ay, the clustering transition. Indeed the configuration space of these models 
gets fractured in a large number of pure states, called clusters in the context of 
random CSP, as explained in Sec. 3.7 and sketched in Fig. 7. The correlation 
decay hypothesis only holds for the Gibbs measure restricted to one pure state, 
not for the complete Gibbs measure. This complication can be handled by the 
cavity method with "replica symmetry breaking" (RSB). It amounts to make 
further self-consistent hypotheses on the organization of these pure states, and 
on the correlated boundary conditions they induce on the tree-like portions of 
the factor graph. Inside each pure state the RS computation holds true, the 
RSB computation is then a study of the statistics of the pure states. Let us 
explain how this is done in practice at the first level of RSB (1RSB cavity 
method). The partition function is written as a sum over the pure states 7, 
Z = Y^-y > where Z 1 is the partition function restricted to the pure state 7 
(recall the decomposition of Eq. (25)). It can be written in the thermodynamic 
limit as Z 1 = e~ N ^, with / 7 the associated (internal) free energy density. 
One further assumes that the number of pure states with a given value / of the 
internal free energy density is, at the leading exponential order, e N ^^', where 
X is called the complexity, or configurational entropy. The latter is assumed 
to be a concave function of /, positive on the interval [/ mm , /max]- In order to 
compute E one introduces a parameter m (called Parisi breaking parameter) 
conjugated to the internal free energy, and the generating function of the Z~ ( as 
Z(m) = Z™. In the thermodynamic limit its dominant behavior is captured 
in the 1RSB potential $(m), 



$>(m) = lim — logZ(m) = inf 

v ; m/3 jv^oo N f 



f-\Mf) 

mp 



(79) 



where the last expression is obtained by a saddle-point evaluation of the sum 
over 7. The complexity function is then accessible via the inverse Legcndrc 
transform of $(m) [235], or in a parametric form 

£(/(m)) = m 2 /3$'(m) , /(m) = $(m) + m$'(m) , (80) 

where f(m) denotes the point where the infimum is reached in Eq. (79). 

The actual computation of $(to) is done as follows. One introduces on 
each edge of the factor graph a distribution Pi^, a (rj) of messages, which is the 
probability over the different pure states that /ii-> a = ?7, where fj-i-> a are the 
messages that appear in Eq. (71). Because Pi_j. a (ry) fluctuates from instance 
to instance, the order parameter now becomes the distribution of Pi^ a (ri) with 
respect to the disorder, which we call V^ 1 '{P) and is solution of a self-consistent 
functional equation written as 

P = G{Py, u ...,Pr, k -u--^Pd,r,-^,Pd,k-i,v,w x ,...,w d ) . (81) 
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Similarly to the RS case the P's are independent copies drawn from pW, and 
the random weights v and w a are also independently generated. The r.h.s. of 
this equation stands for: 



Z{{Pa,i}' V A w a},m) 

J [] dPaAVa4)S(V-9({V*M<{Va,i},V,{Wa}) m , (82) 



ae[l,d] 
»e[i,fc-i] 



with g and z defined above in Eq. (76), and to is the Parisi parameter. From 
the solution of this equation one computes the 1RSB potential $(m) via an 
expression similar to (77), namely 



-m/3$(m) 



E 



log 



J Yi dP a j(r] ati ) z v ({r) a .i},v, {w a }) T 



\ 



a{k - 1)E 



ae[l,d] 
ie[i,fc-i] 



log 



dPi(r)i)zf({r)i},w) 

ie[l,k] 



(83) 



with the functions z v and Zf defined in Eqs. (78), (79). As in the RS case this 
expression is variational, the implicit dependence of P' 1 ) on the various param- 
eters (to in particular) can be discarded when computing derivatives. Various 
physical situations translate in different behaviors of the 1RSB equations. 

• It can happen that only trivial solutions of (81) exist, i.e. the P's are 
supported on a single value of rj. It is then easy to check that this random 
T] obeys precisely the RS equation (75), and that the 1RSB potential <I>(to) 
is equal to the RS prediction for the free energy of Eq. (77), for all values 
of to. This is the translation of the existence of a single pure state, in 
which case the whole 1RSB machinery reduces to the RS case. This case 
is realized at high temperatures/low connectivity parameter a, i.e. on the 
left of the line Td(a) of Fig. 8. In more physical terms it corresponds to 
a "liquid" phase, for instance the high temperature phase T > Td of the 
fully connected p-spin model described in Sec. 3.4. 

• If on the contrary non-trivial solutions of the 1RSB equations appear, 
one has to investigate them more carefully in order to obtain the 1RSB 
prediction for the free energy density. From the definition of Z(m) one 
would naturally take $(1) for it. This is indeed the case, provided the 
corresponding complexity S(/(m = 1)) is positive, in other words if 
/(to = 1) G [/min, /max]- The physical interpretation is that an exponen- 
tial number e N ^ f ^ of pure states contribute to the Gibbs-Boltzmann 
measure, each with an internal free energy density f(m = 1). It turns 
out that in this case the prediction $(m = 1) coincides with the RS one; 
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this form of replica symmetry breaking is not seen in the thermodynamic 
properties of the model, yet it has drastic consequences on its dynam- 
ics [236, 179]. Such a phase is usually called for this reason a dynamic 
1RSB (dlRSB) phase, and is realized in the fully connected p-spin model 
in the intermediate temperature regime T c < T < Td, or more generi- 
cally for diluted mean field spin glasses in the part of their phase diagram 
enclosed by the lines Td(a) and T c (a) (see Fig. 8). 

• If there are non-trivial solutions of the 1RSB equations, but with E(/(l)) < 
0, one has to find the value m s <G [0, 1) which solves the equation S(/(m s )) = 
0, i.e. f(m s ) = ,/min, and the 1RSB cavity method predicts that the free 
energy density is precisely this value / m j n . In more physical terms this is 
called a condensed (or true 1RSB) phase, many pure states exist in the 
system yet only the sub-exponentially numerous ones with / = /,„;„ do 
contribute to the thermodynamic behavior of the system, as for instance 
in the T < T c phase of the fully connected p-spin model. 

As far as the free energy density is concerned, one can give a general formula 
that encompasses and summarizes the three cases above: it is given by the 
maximization of the potential $(to) with respect to m, 

f = — 4 lim — EflogZl = max <6(m) . (84) 
J Pn^ooN 1 8 J me[o,i] 17 v ' 

5.1.5. Population dynamics 

In general there is no hope to find an analytical solution of the cavity equa- 
tions, neither at the RS level (75) nor at the 1RSB level (81). These equations 
are however amenable to a numerical resolution in a relatively simple way. 

Let us first consider the RS case. The variables rj are probability distributions 
over the discrete space X, each of them can thus be represented by \X\ — 1 real 
numbers (thanks to their normalization). In particular for Ising spins a single 
real is enough, that can be interpreted as an effective magnetic field acting on 
a spin. In the RS equation (75) the unknown is V(rj), a probability distribution 
over such effective fields. A convenient way to represent it numerically [237, 152] 
is to use a sample, also called population, of representative fields, i.e. to write 

-, JVext 

V{ V) = ttE%- fh) ■ (85) 

1=1 

This representation thus uses a number A4xt of rcprcscntants rji, each of them 
encoded as a single real for Ising spins, or \X\ — 1 real numbers in the generic case, 
and is obviously more and more accurate as A4xt gets larger. As the equation 
(75) has the form of a fixed point condition, it can be solved by iteration: 
one starts from an arbitrary initial assignments of the sample {r^}, plugs the 
representation (85) in the r.h.s. of (75), and constructs a new set {ij^} that 
represents the l.h.s.. This is simply done by repeating Af cx t times the following 
steps: 
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• draw an integer d from the law qd. 

• draw d{k — 1) integers . . . ,jd,k-i uniformly at random in [l,A/" C xt]- 

• draw a random vertex weight v and d random interaction weights w\, . . . ,Wd- 

• compute rf from Eq. (76), where rj a ,i is taken to be the j a ,i'th element of 
the current population representation of V . 

The 7V ex t elements rf thus generated form a sampled representation of the l.h.s. 
of (75); this process can be iterated, i.e. this new representation can be plugged 
in the r.h.s., and so on and so forth until convergence towards the fixed point of 
(75) is achieved. Then the expectation values over T 3 , as for instance in Eq. (77), 
can be simply computed by taking an empirical average over the sample of the 
jV ex t representants of V, which gives access to all the physical obscrvables of the 
model, at the replica symmetric level. 

Let us now discuss the generalization of this method to the 1RSB level. The 
equation (81) has exactly the same structure as (75), and one can thus follow 
the same strategy as above, with a population of A4xt representants Pi of the 
distribution The last step in the algorithm explained above becomes 

• compute P' from Eq. (82), where P a ^ is taken to be the j a ,i'th element of 
the current population representation of 

This step itself deserves some explanation. The additional difficulty is that Pi is 
itself a distribution over fields; but this can be handled in a similar way, by rep- 
resenting each Pi by a sample of Mnt fields rji^, with i e [l,A/" ex t], i' G [l,A/lnt]- 
Now to generate a representation of P' from Eq. (82) one has to construct 
A/lnt fields r], by extracting the r) a ,i from their respective distributions P 0) , and 
computing ?/ = g({r] a ^}). However these A/l n t fields should not be given an 
equal importance in the representation P'\ they have to be weighted according 
to the replica symmetry breaking weighting factor z{{r] a ^} 1 v 1 {w a }) m - Several 
reweighting schemes can be used to perform this task [152], and this is a well- 
studied problem in the field of statistics where this kind of representation is 
called "particle approximation" [238]. A simple idea to perform this reweight- 
ing will be given for a similar equation in Sec. 5.2.2, we refer the reader to the 
original literature for more details on this issue. 

To summarize, the generic RS cavity equation is handled numerically by 
a population (of J\f cx t elements) of fields (each corresponding to |x| — 1 reals, 
i.e. a single one for Ising spins), while the resolution of the 1RSB equation 
involves a population of A/" ex t populations of A/lnt fields. Higher levels of replica 
symmetry breaking [32] can be formally treated by increasing the number of 
generation in this hierarchical construction, but in general the numerical cost 
of their resolution increases too fast to allow one to go beyond the first step. 

There are however some cases, to be encountered in the following, where 
one level of complexity disappears. For models defined on regular random (hy- 
per)graphs, the distribution qd is concentrated on a single integer; then, if the 
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vertex and interaction weights v and w are not random or arc sufficiently sym- 
metric, the RS equation (75) (resp. the 1RSB equation (81)) admits a solution 
where V (resp. V^) is concentrated on a single field rj (resp. on a single distri- 
bution P). In other words the external size A/" ex t can be reduced to 1 in these 
cases, which greatly reduces the numerical cost of the resolution of the equations 
(in particular in the quantum will be discussed below); this situation is 

often termed factorized in the mean field spin glass literature. 

5.1.6. Analyzing a mean field spin glass model 

The phase diagram of mean field spin glass models, as a function of the 
temperature, magnetic field, connectivity, or other control parameters, can vary 
qualitatively from model to model. We will not attempt here to provide a 
general classification, but only propose a flowchart that one should follow when 
confronted with a new model, for each point of its phase diagram. As a first 
step one should solve the RS equation of the model, and compute the physical 
observables associated to it. Then the validity of the RS hypothesis should be 
tested; in some cases its violation is apparent from the inconsistency of the RS 
results (a negative entropy for instance), but not always. In any case the 1RSB 
equation has to be solved, first with the breaking parameter m set to 1. Two 
cases can then appear: 

• if there are no non-trivial solutions of the m = 1 1RSB equations, the RS 
results should be conjectured to be valid, the system is in a liquid phase. 

• if there is a non-trivial solution of the m = 1 1RSB equations, then one 
has to check the sign of the associated complexity. 

— if S(/(m = 1)) > 0, the system is in a dynamic 1RSB phase. The 
RS prediction for the free energy is correct, yet the Gibbs-measure is 
split on an exponentially large number of pure states (clusters) , and 
the dynamics is non-ergodic. 

— if £(/(m = 1)) < 0, the system is in a true 1RSB glass phase, with 
a sub-exponential number of pure states dominating the equilibrium 
measure. It is then necessary to make a study as a function of m £ 
[0, 1], and to find the value m s where the complexity vanishes. The 
1RSB prediction for the free energy is then /(m s ). 

As a matter of fact this analysis can be further complicated by the co- 
existence of multiple solutions to the RS or 1RSB equations (that should be 
discriminated by comparing the free energies they yield), and by the instability 
of the 1RSB solutions towards higher levels of replica symmetry breaking. This 
second point is particularly difficult to handle in the case of diluted models [239]; 
however, in fully connected models where any level of RSB can be solved, it is 
found that the 1RSB results are quantitatively very good approximations to 
exact ones. 
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5.2. The path integral quantum cavity method 

5.2.1. Path integral representation of discrete quantum models 

We shall now define the quantum version of the models studied in the follow- 
ing. To the space X N of classical configurations and energy function E(a) we 
associate the Hilbert space spanned by the vectors {|cr)|<7 e X N }, and an oper- 
ator Hp (we keep here the notations of Sec. 4), diagonal in this basis, according 
to Hp = J2 a E(&)\cr)(a\. We then define the quantum Hamiltonian 

N 

H = H P -J2^ l T l , (86) 

i=l 

where in the second term the operator acts on the i-th variable only, according 
to 

(a\fi\a') =T au< n<W;. • (87) 

Without loss of generality we assume that the matrix T of order j^l has van- 
ishing diagonal elements (these can be incorporated in the classical part Hp). 
We make the further hypotheses that all its off-diagonal elements are real non- 
negative (this is crucial to avoid the "sign problem" ) , and that T is symmetric 
(this ensures the Hermitian character of H). The simplest example is provided 
by the spin 1/2 case, with X = {— 1,+1}; the above basis is taken to be the 
eigenvectors of the Pauli matrices af , and one can take Ti = af , the parameters 
T,; are then the transverse fields of the model. The quantum statistical mechan- 
ics study of the model amounts to the computation of the partition function 
Z = Tre _,3fl . The additional technical difficulty, with respect to the classical 
models, arises from the non-commutativity of the two terms in H. The standard 
way to handle this difficulty is to introduce a path integral representation of the 
partition function, of the form: 



~[Der. l v l ((T l ) 

i=i 



2l(0)=<x(/3) 

Let us precise the notations we introduced. Here and in the following bold 
symbols will be used for functions of an "imaginary" time t £ [0, j3]; in particular 
here <Ji is a piecewise constant function <7j(t) : [0,f3] — > X. The integration 
measure Deri is decomposed as a sum over the number n of discontinuities of 
the function <Ji(t), the times t\ < ■ ■ ■ < t n at which they occur, and the values 
of, . . . , <t" the function takes on the n + 1 intervals [0, ti], [t±, £2], ■ • ■ , [t n , 0\- 



00 



DtTi = Y, / dti / d * 2 • • • / d< " ■ ( 89 ) 



n=0 CT o. 



Jti 



For such a function the weight Uj(o"j) reads 



M = (Pi) n U T ^tJ)Mtt) = ^ n U T ai-\ai S ( 9 °) 

j=i 1 3 i=i 
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as the diagonal elements of T are supposed to vanish the only contributing paths 
are those with a\~ x 7^ o{ . Such a path integral representation can be devised for 
any matrix element of e _/3 ^, by simply fixing the initial and final configurations 
cr(0) and g_({3). Here we let them free, under the condition <r(0) = u_{0), to 
compute the trace of e~@ H . 

There are several ways to obtain such a path integral representation. The 
most pedestrian one is to use the Suzuki- Trotter formula, decomposing the 
imaginary time interval [0,/?] in N s slices, 

e Xi+x 2 = Hm / ^-x le -fcx 2 \ /gjN 

N s ^oo V / 

with the two non-commuting operators X\ = —/3Hp and X 2 = /S^-TiTj. One 
can then introduce N s representations of the identity between the terms of the 
product, as sums over the configurations a(t = a/3/N s ) with a a discrete time 
index. The continuous time limit N s — » 00 then yields the path integral (88). 
A maybe more elegant way to obtain this result is to use the following operator 
identity, 

e x 1+ x 2 dt l dh / dtp e HX, £ 2 e (fa-ti)*i x 2 ...X 2 e^-Wi , 

with the same values of as above; one then inserts p representations of the 
identity in the p-th term of the sum, rescales the time integrals and reorders the 
summation according to the number of times each flipping operator Tj is picked 
in the expansion of X 2 . Finally these path integral representations can also 
be handled in a mathematically rigorous way, see for instance [240, 241, 242, 
243, 244, 245, 246], by starting from Poisson point processes for the candidate 
times of discontinuities of the variable trajectories, a distribution which is then 
properly biased by the classical energy terms. 

5.2.2. Representation of the cavity messages 

The path integral representation of the partition function given in Eq. (88) 
is valid for any classical energy function E{q). Suppose now that E is decom- 
posed as a sum of local interaction terms, according to Eq. (69). The quantum 
partition function (88) can then be rewritten as 

M N 

z = E ^(0)^0?) n^-C^n ' ™a(<L 9a ) = e- f ° dt£a( ^ m - 
&ex N a=1 i=1 

(92) 

where we introduced X the space of periodic piccewise constant functions from 
[0, 0] to X ', and used the notation X)<t eA 5 as a s y non y m of the integration J Deri 
defined in Eq. (89). This notation was chosen to emphasize the similarity with 
the classical partition function (70): the quantum computation is reduced to a 
classical one, the cost to be paid being the replacement from a discrete variable 



<7j in X to a function (trajectory) er; = {<7i(t)\t e [0,/3]} in X as basic degrees 
of freedom. The partition function can also be interpreted as the normalization 
constant of a probability measure over X N , namely 

j M N 

Mq(ct) = II WaioLdjJl^i) ■ ( 93 ) 

Note that apart from the change of the nature of the degrees of freedom, the 
"spatial" structure of the interactions w a encoded in the factor graph is the same 
in the classical and in the quantum case. In particular as soon as the classical 
energy part of the quantum Hamiltonian falls into the category of models that 
can be solved by the cavity method (i.e. sparse random graphs), then this is 
also true for the quantum problem. This observation was first exploited in [223] 
to study the spin 1/2 spin glass on the Bethe lattice, within a finite number 
of Suzuki- Trotter slices. The formulation of the quantum cavity method in 
continuous imaginary time was then presented in [224] at the RS level, for a 
ferromagnetic model, and [232] for the Bose- Hubbard model of bosonic particles. 
Results at the 1RSB level were given in [48], and a complete exposition can be 
found in [198] for interacting particle models. 

The general structure of the quantum cavity method is thus exactly the same 
as the classical one exposed in Sec. 5.1. In the rest of this section we explain 
the additional technical points that arise when going from X to X as the base 
space for degrees of freedom. First of all one has to find an efficient way to 
represent the probability distributions r](er) over X, which are the basic objects 
of the method. As should be clear from the discussion of Sec. 5.1.5, the simplest 
way to do that is to approximate them by a weighted sample representation of 
a large number A/traj of elements of X , namely 

JVtraj A/traj 

n{a) = aU) S(tr - <r w >) , with ]T a® = 1 . (94) 
i=i 3=1 

Each representative trajectory rr^' is itself encoded in a compact way by the 
number of discontinuities it contains, their times of occurrence, and the constant 
values of the function between the discontinuities. Sampling an element from rj 
means drawing a number j G [1 , A/traj] with probability a"', and returning the 
trajectory cr^\ 

Secondly one must devise a way to implement the quantum equivalent of 
Eq. (76), 

z({Va,i}, T, {e a }) 

e (n n^ j ° dtEaWt)iff " i(i) . (95) 

f^.KS'V \ a ' 1 ' a=1 
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with the quantum variable weight defined in Eq. (90), 

n 

B ( ff ) = rn II T ^) > ( 96 ) 

3=1 

n denoting the number of discontinuities of <r, at times t% < ■ ■ ■ < t n . To 
rewrite this equation under a more convenient form we shall introduce time- 
varying "longitudinal fields" h, which are functions h a {t) of an index a G X 
and of the imaginary time t. Indeed the dependency on <r of the energy terms 
in Eq. (95) can be conveniently expressed in terms of such a field defined by 

d 

h{{(T a ^},{e a }) : h a (t) = - e a (a, a a ,i(t), a a ,k-i{t)) ■ (97) 

a=l 

We can then rewrite (95) as 

p(«r|I\ £({*„,<}, { e .})) Z( y ({ ^ } ; {£ ° })) , (98) 

where we defined 

p(a-|r,h) = 1 —v(<r)eS° ' dth -wW , (99) 

z(r,h) 

In this way p(er|r, ft.) is a well-normalized probability distribution over X . Sup- 
pose that one is able to sample from this distribution, and to compute the 
associated normalization Z(T, h) (we shall see in a short while that this is in- 
deed possible). Then the resolution of Eq. (98) is at hand. This amounts in 
fact to the generation of a sampled representation of its l.h.s., assuming the 
knowledge of the r.h.s., in particular the ability to draw the trajectories from 
the probability laws rj a i on X. To construct the representation of the l.h.s., one 
has to repeat At ra j times independently, for j £ [l,A/t ra j], the following steps: 

• draw the {(T 0| i}Sj ^ ^ from their respective distributions rj a a- 

• compute the field h according to Eq. (97). 

• extract a trajectory cr^' from the law p(-|r, h) and set 0"' = Z(T, h). 
Once these steps have been performed A/t ra j times, normalize the new weights, 

« 0) <" -m nr-^ ■ (100) 

a (l) + . . . + ffl (Mraj) V > 

A moment of thought reveals that this is indeed the correct algorithmic trans- 
lation of Eq. (98). 
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5.2.3. Path generation 

We are thus left with the problem of generating a path according to the law 
p(cr\r, h) defined in Eq. (99) and of computing the normalizing factor Z(T, h). 
As all trajectories cr ai are piecewise-constant, this will also be the case of the 
relevant realizations of the field h. Let us call p the number of discontinuities 
on [0,/3] of h, that occur at times < v- 1 ' < ■■■ < v?) < /?, and denote 

the values of h(t) on the time intervals [0,*W], [tW,tM], 
[t&\0]. We also denote A< ) = tW, \W = - i«, . . . , A (p) = 0- t<» the 
duration of these intervals. Let us introduce some further notations; we consider 
a \X\ dimensional Hilbert space spanned by {|er)|cr € X}, on which we define 
an operator H(T, h) by its matrix elements, 

(a\H(T,h)\a') = -h a 5 a , a , -TT a . a , . (101) 

We shall write W(T, h, A) = e~ XH< - r ' h ^ its associated propagator on an interval of 
imaginary time of length A, and W(T, h, X) at<7 > = (a\W(T,h,X)\a') the matrix 
elements of the propagator. It is then possible to prove that the sought-for 
normalizing factor Z(T, h) reads 



Z(T,h) = Tr 



]iy(r,ftW,A (,) ) 



(102) 



This is a computationally affordable expression: it requires diagonalizing p ma- 
trices of (small) dimension \X\, exponentiating them and multiplying them to- 
gether. Finally the process of generation of cr with the law p(-\T,h) can be 
implemented as follows: 

• draw the values , . . . , that cr assumes at times 0, v- 1 ' , . . . , t^ p ' . 

• on each of the p+1 intervals [v^ , t^ l+1 '], draw a trajectory representative of 
the evolution W(T, h^,X^) in a constant field hfr' , with boundary condi- 
tions <r(t«) = o-t l \ ( j(t( l + 1 )) = (we set t<°) = and = a^). 

More precisely, the first step consists in extracting these p+1 values from the 
joint law 

p(cr(°), . . . , >>) = —L^W(T, h<-°\ A(°)) CT (o, iCT(1) W(T, fcW, X^Uv^ ■ ■ ■ 
Z{T,h) 

...W{T,h^\X^) a{p) ^ a) . (103) 

This can be easily done by first drawing a^ 1 from its marginal probability, then 
cr*- 1 ) conditioned on the value of and so on until has been extracted. 

The procedure to follow for the second step is more apparent once an integral 
equation on W is written: 

W(T, h, X) a , a > = e Xh "8 a , al +T [ dt e tK V T a ^W(T, h,X- t) a ,^ a , . (104) 

Jo 
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In terms of the path integral representation of W, the two terms in this equation 
represent respectively the contribution of a constant path (possible only if the 
boundary conditions are the same at time t = and t = A) and of a path 
whose first discontinuity occurs at time i, where o~(t) jumps from a to a" . In 
consequence, the procedure to draw a path from o~(t = 0) = a to a(t = A) = a' 
in presence of a constant field h reads 

• if a = a' , with probability e xha /W(T, h, X) a . a , exit with the constant path 
a(t) = a Vt 6 [0,A]. 

• otherwise 

— draw a random time u G [0, A] with the cumulative distribution 

; A at e*^ £ CT „ r CT>CT "W(r, h,\ — t) a n^, 



— draw an element a" with probability 

T (Ti<r «W(r,/i,A-u) <T « (7 . 



E^T a>(T ,„W(T,h,X- 



(106) 



— set cr(t) = a for t € [0, u], and call recursively the same procedure 
to generate the path on [it, A], with boundary conditions <j(u) = a", 
a(X) = a'. 

This procedure for the generation of imaginary time paths in presence of a 
constant transverse operator T and piecewise constant longitudinal fields was 
presented in the case of spins 1/2 in [224]; its recursive nature has the advantage 
of making it a rejection-free method, the inconvenience being the necessity to 
draw random variables from rather complicated distributions (105) for the in- 
terval of times between spin flips. An alternative method was proposed in [41]: 
the time intervals between spin flips are drawn from exponential distributions 
with well chosen averages, which is much easier, but the price to be paid is a 
rejection if the parity of the number of spin flips on the interval [0, 0] does not 
satisfy the boundary conditions on <r(0) and er(/3). The rejection rate can in 
particular become quite high if er(0) ^ for low transverse fields. The two 
methods can actually be combined to gain both their advantages, by drawing 
for anti-periodic trajectories the time of first flip from the recursive method, 
then continue with the rejection one. 



5.2.4- Discussion 

The Path Integral Quantum Cavity (PIQC) method described above is an 
exact way of dealing with quantum spin models on sparse random graphs. The 
analysis of such models proceeds along the same lines as explained in the classical 
case in Sec. 5.1.6, i.e. via an interpretation of the solutions of the RS and 1RSB 
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equations. In Sec. 6.2 and 6.3 we shall present explicit results on two models of 
random constraint satisfaction problems obtained in this way. 

It is however important to mention the limitations of the method. As already 
explained at the beginning of this section, it can only handle quantum systems 
that do not suffer from the sign problem, and it is a hnitc temperature method; 
the ground state properties are necessarily obtained by an extrapolation to zero 
temperature. Moreover the numerical resolution of the quantum cavity equa- 
tions is a numerically costly task. The generic quantum 1RSB case involves 
indeed a representation by A/" e xt x Mnt x A/t ra j imaginary time trajectories 4 
(recall the discussion of Sec. 5.1.5 in the classical case and the additional pop- 
ulation level due to the quantum nature of the model explained in Eq. (94)). 
Fortunately for factorized models (with regular degrees) this is reduced with 
■A/"ext = 1, see the end of Sec. 5.1.5. In any case the memory available on present 
days computers limit the numbers Mnt and At ra j to relatively small values (ex- 
amples will be given on concrete cases in Sec. 6.2 and 6.3). This induces both 
systematic deviations of the empirical mean from the exact value and noise in 
its estimation; extrapolations to ./Vint, A/traj oo via finite size analysis can 
however be performed to reduce these effects. A specific difficulty comes from 
the weighted representations of probability distributions used in Eq. (94) for 
instance; one must take care by resampling methods of the tendency that these 
weights have to flow towards very inhomogeneous repartitions, which leads to 
situations where the number of effective representants of the distribution be- 
comes much smaller than A/traj [238]. 

In the following sections we shall describe alternative approaches to quantum 
models on sparse random graphs that, even if approximate, allow to bypass some 
of these limitations. 

5.3. Operator quantum cavity methods 

In this section we shall describe alternative formulations of the quantum 
cavity method that do not make use of the path integral formulation but work 
directly with quantum operators [247, 225, 226, 227, 228, 229, 248, 230]. These 
approaches have been sometimes called "quantum belief propagation" but we 
will refer to them here as Operator Quantum Cavity (OQC) methods. They 
all share common features and ideas whose connections are still only partially 
understood. They represent approximated methods, and their level of accuracy 
is not completely controlled yet. However, compared to path integral methods 
they have the important advantage that the T = limit can be taken explicitly. 
Moreover, the cavity messages are represented as finite matrices, therefore there 
are no sampling errors, unlike in the PIQC where the messages are represented 



4 An alternative approach, that will not be further discussed here, consists in a systematic 
perturbative expansion in the transverse field T; any finite order of the expansion can be 
expressed in terms of the classical cavity computation, thus strongly reducing the numerical 
cost with respect to the fully quantum approach. This however does not give access to non- 
pcrturbativc effects like phase transitions. 
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through finite samples of probability distributions over an infinite-dimensional 
space. 

For the sake of simplicity, we will present these methods in the simpler case 
of an Hamiltonian that is the sum of two-body interactions: 



H 



(107) 



The sum over runs on the links of a regular lattice of degree c. We 

will mainly focus on the case c = 2 of a one dimensional chain, and c = 3 
with the underlying lattice being a 3-regular random graph. Moreover we 
will consider the case of an Ising model in a transverse field for which Hi j = 
—Jijafa? — (Ti&f + F ja?) / c. The generalization to more complex Hamiltonians 
is straightforward. 

5.3.1. Operator cavity messages 

We start our presentation by following the derivation of [227] and consid- 
ering for simplicity a finite one-dimensional chain with open boundaries. The 
quantum partition function is 



Tre-^=Tr 



Qe 



-13H 2 



© 



e 



(108) 



where e A e B = e A+B . As in the classical case, we can define an operatorial 
message that acts on the Hilbcrt space of spin i only: 

1 



"Tri, 



i-i e 



~P 5Ifc=i Hk,k+1 



(109) 



where the normalization is determined by Tr $ ?7i-n+i = 1. 

We can derive an approximate recurrence equation for these messages by 
following the same steps as in the classical case: 



J7i_M + i oc Tr | 
= Tr. i _ 1 {Tr 1 



-PT. 



e 



k=l H k,k+1 q e 



-pHi- 



Tr 



Tri ... i_ 2 e 



-PYXJi Hk,k+i 



© e" 



(110) 



cx Tri_i [Tfi-x-^% © e" 



-PHi 



and the proportionality constant is determined by normalization as in the clas- 
sical case. The crucial point is that, unlike in the classical case, here we made 
an approximation when we changed the position of the square brackets moving 
from the second to the third line of the above equation. 

Indeed, consider a system made of three parts a, b, c and operators H a b, 
H\y C , acting only on a® b and b®c respectively. Due to quantum entanglement 



Tr, 



e -PH a , b Q e -fm b 



Tr„e" 



-PH a 



© e 



-PH b 



(111) 
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However, the above equation is an equality if the "conditional mutual informa- 
tion" I(a : c\b) = S(a, c) + 5(6, c) — 5(6) — S(a, 6, c) vanishes (here 5 is the von 
Neumann entropy), indicating that all correlations between a and c are medi- 
ated through 6 (as in the classical case). It has been argued that this condition 
holds when the region 6 is sufficiently "thick" [227]. The problem is that in 
Eq. (110) the region 6 coincides with a single spin, 6 = {i — 1}. 

This observation motivates the introduction of new messages, that are oper- 
ators on the space of spins {i — £ + 1, • • • , i}. Repeating the above derivations: 



1 



oc Tr ;_ f { Tr i, 



i-l-l 



0e 



Tr, 



{[Tri,-,* 



-ie 



~P 12k=l Hk,k 



©e 



(112) 



The crucial difference is that now the region 6 = {i — £,..., i — 1} has thickness 
t and one can hope that the error is much smaller. An argument in favor of this 
has been discussed in [227] . The drawback is of course that now the messages 
are operators acting on I spins, and therefore they have to be represented by 
matrices of size 2 e . 

The generalization of this procedure to a tree is straightforward. Let us call 
Ti-tj the partial tree rooted at i obtained by cutting the link and d(i,j) 

the distance on the tree between i and j. The message from i to j is defined as 

1 



i->3 



(113) 



and we get as in the classical case: 



'/ 



oc Tr 




kedi\j 



(114) 



Here, the messages are operators acting on 1 + (c — 1) + (c— 1) 2 + ■ ■ • + (c— l) e 1 
spins, so they must be represented by matrices whose size =o( c !) grows 
much faster than in the one dimensional case. 

With similar reasonings one can obtain the approximate expression for the 
free energy, which is exactly the same as in the classical case (here specialized 
to a system with two-body interactions only) , with sums replaced by traces and 
the normal product replaced by the product: 



-pF 



E 



log * - ^2 lo S - 



■13 i 



za = Tr 



(115) 



rp I (i) „ (I) ^ - 



13 Hi 
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5.3.2. Explicit equations for single- spin messages 

Let us now consider more explicitly the above OQC formulation on a tree 
with £ = 1. In this case the messages are operators on a single spin, i.e. 2x2 
Hermitian matrices normalized to have trace 1. We can parametrize them by 
two local fields: 



1 



(116) 



omitting a term proportional to erf that vanishes by symmetry. Equivalently 
we can describe the message f)i-^j in terms of the magnetizations 



Tr i (o : fr/ i ^. i ) 



: tanh 



b 2 ■ 



: tanh 



b 2 ■ 



Plugging this in Eq. (114) with £ = 1 we obtain 
which can be recast in the following form: 



Tr^ gaA3 (afe-^) 



and similarly for , where 



kGdi\j 
k£di\j 



(117) 



(118) 



(119) 



(120) 



is an effective Hamiltonian acting on spin i and its neighbors (except j). It- 
eration of these equations then requires at each step the diagonalization of a 
Hamiltonian acting on c spins. Note that taking the T = limit is straightfor- 
ward and simplifies the computation, because in this case we only need to find 
the ground state of if e ff- 

One can actually take a different approach and substitute Eq. (116) in the 
free energy Eq. (115), obtaining then a function of the set of all fields b^j and 
hi-yj. One can then derive equations for these fields by imposing stationarity 
of the free energy with respect to variations of any field, as in the classical 
case. However, because the OQC is only approximate, the stationarity equations 
do not coincide with the equations obtained from cavity iteration, Eqs. (117), 
(119), (120). It can be shown on specific examples (e.g. the ferromagnetic case 
Jij = J and Ti = V) that imposing stationarity of the free energy is slightly 
more accurate than the iteration scheme. 
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Let us also mention a further approximation that has been proposed in [229, 
248, 230], which amounts to replace the operators erf,?! in Eq. (120) by their 
averages m^^mj^ in Eq. (117). One thus obtains the following equations: 



c — 1 

Oi— >j — 1 i 



kedi\j kedi\j V "fc^i + L 



fc— >i 



These are closed and relatively simple equations for the fields h^j and have been 
exploited in [229, 248] to obtain detailed information on a disordered system 
that would have been extremely hard to obtain from the numerical solution 
of the OQC or PIQC equations. Additionally, it is clear from these equations 
that one can take the /3 — > oo limit without problems just by dropping the 
hyperbolic tangent term. A drawback of this approach is that these equations 
are approximate, even in the classical case Ti = 0. It has been argued in [229, 
248] that they become exact for c — > oo, see [230, 249] for a detailed discussion 
of this delicate point. 

5.3.3. Relation with the PIQC 

The OQC has been introduced in [229, 248], independently from [227], to 
study the metal-insulator transition in disordered superconductor and later used 
in [230] to discuss the properties of disordered fcrromagncts. The derivation 
of [229, 248, 230] starts from the PIQC formulation and makes a simple ansatz 
on the functional form of the distribution of imaginary time trajectories. In 
turn, this can be reinterpreted as an ansatz over the Hamiltonian governing a 
reduced part of the system, consisting of neighboring spins, and gives back the 
OQC. 

The PIQC leads to the following equation (which is the specialization of the 
treatment of Sec. 5.2 to Ising spins, see also [224]): 

7? i _v J -(<7 i ) = — II f Dtr k ri k -n(tr k )e J «ti' r *V'"'W dt , (121) 
Zl ^ 3 kedi\j J 

where |<Tj| is the number of spin flips in the imaginary time trajectory cr^. In 
order to simplify the solution of these self-consistent equations, in [229, 248, 230] 
it was suggested to consider the following ansatz: 

rii^iiai) oc (^ J .) kil e / " Wi(t)d< ■ (122) 

Once inserted in the right hand side of Eq. (121) this ansatz doesn't give back 
in the left hand side a message of the same form. However one can take its 
"projection" over the distributions of trajectories described by (122), by fixing 
the new fields and 6j->j in such a way that the expectation values of af 

and erf on the two sides of Eq. (121) are the same. Not surprisingly, it is easy 
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to show that this procedure gives back 5 the same equations as the OQC for 
1 = 1, Eqs. (117), (119), (120). It was shown in [230] that this approximation 
gives quite good results when compared with the exact PIQC solution, and the 
quality of the approximation increases with increasing connectivity c. 

For I > 1, the connection between OQC and PIQC is less obvious. We will 
not discuss it in detail, but roughly speaking the idea is the following. The 
Markovian ansatz in Eq. (122) neglects all imaginary time correlations in the 
path integral description of spin i. Therefore, a more refined ansatz would 
include, for instance, a Gaussian term J Q dtdt'Gi^j ;{t — t')ai(t)ai(t') in the 
exponent [223]. In presence of such a term, the PIQC equations cannot be 
cast in an operator formulation using only local operators. The reason is that 
these imaginary time correlations are obtained by tracing out the neighboring 
spins. In the PIQC representation, this could be represented by considering a 
Markovian ansatz acting not only on i but also on a neighboring shell of size t, 
and then integrating out the neighbors to obtain an imaginary time correlated 
message on spin i. In the OQC language, this should correspond indeed to an 
operator message acting on spin i and a set of neighbors. We conclude that 
messages with £ > 1 in the OQC should roughly correspond to adding some 
imaginary time correlations in the PIQC. This is very reminiscent of what is 
done in dynamical mean field theory where imaginary time correlations are often 
represented by an Hamiltonian thermal bath of phonons [250] . 

5.3.4- Discussion 

OQC [247, 225, 226, 227, 228, 229, 248, 230] (or Quantum Belief Propaga- 
tion) is a very promising approach to the solution of spin glass models on locally 
tree-like graphs. First of all, this method is not affected by the "sign problem" 
and therefore can be applied to Hamiltonians that do not admit a path integral 
representation with positive weights (e.g. the QSAT problem [55]). Another 
important advantage is that for a given £ the cavity messages are finite matrices 
that can be parametrized by a finite set of real numbers. The accuracy of this 
representation is only limited by machine precision, unlike in the case of PIQC 
where sampling introduces systematic numerical errors and noise. For a given 
£, the limit T = can be taken easily by replacing everywhere the traces at 
finite temperature by a ground state average. 

Its main drawback is that it is an approximate method: its accuracy is 
expected to increase with the size of the block £. If one requires a given accuracy, 
then the block size must be increased when decreasing T and £ — > oo for T — >• 
[227] (however, there is some hope to combine OQC with local renormalization 
group methods to avoid this problem [227]). At the same time, for a fixed block 
size, the limit T — > exists, is simpler to handle than the finite T computation, 



5 Actually, there is a slight difference due to the fact that in the OQC formulation above 
we chose to symmetrize the local Hamiltonian Hij. The PIQC leads naturally to a non- 
symmetric formulation where Hij = —Ji,j&i&j — Fj<Tj/(c — 1). This difference should not 
be crucial, especially for large c where approximation (122) is better justified. 
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and should provide qualitatively correct results [229, 248, 230]. 

The other important problem is that it requires the numerical diagonaliza- 
tion of matrices of large size (especially on a tree with large c, or if the local 
variables are not Ising spins, or if there are many-body interactions). This is 
for the moment a strong limitation to its applicability to interesting problems 
such as XORSAT or the coloring problem as we will see in Sec. 6. 

To conclude this discussion, it is important to stress that OQC can lead to 
rigorous bounds on the true free energy of a given problem. Indeed, in [228] 
it was shown, based on the strong subadditivity property of the von Neumann 
entropy, that slightly modified OQC equations can lead to a lower bound of 
the free energy. Similar results were first derived in the context of classical 
systems [251]. For the moment, this variational technique has been only applied 
to low dimensional systems [228] . Its generalization to random graphs with tree- 
like geometries seems a challenging problem. 

5.4- Variational quantum cavity methods 

A promising approach to overcome the sign problem and to investigate di- 
rectly the T = limit consists in using the cavity method to optimize vari- 
ational wavcf unctions. We will refer to this approach as Variational Cavity 
Method (VQC). These methods are based on the well known fact that, given 
a system described by a Hamiltonian H, the expectation value of the energy 
over an arbitrarily chosen wavefunction \ip v ) provides an upper bound for the 
true ground state energy. We present in this section two different attempts in 
this direction. For simplicity, as in Sec. 5.3 we will focus on the transverse field 
Ising Hamiltonian, H = J2(i j\ Hi-j with H{ t j = —Jijafaj — iVidf +Tj<jj)/c on 
unidimensional chains (c = 2) and 3-regular graphs (c = 3). 

5.4-1- Optimization of Jastrow wavefunctions 

In [231] a very simple trial wavefunction was considered: 

(a\i/i v ) = -^ei b ^ + i ^<-> K W* , (123) 
V Z 

where the constant Z is determined by normalization and bi, Kij are real num- 
bers (it is shown in [231] that adding an imaginary part only increases the 
energy) . Then, the square of the wavefunction is the Gibbs probability measure 
of a classical "auxiliary" system, that turns out to be a classical Ising model 
with couplings and random fields bi at unit temperature. Such a varia- 
tional wavefunction (often called Jastrow wavefunction) has been widely used 
in the study of quantum systems and for many problems it works as a very good 
approximation [252]. 

Once that the wavefunction is chosen one needs to express the expectation 
value of H as a function of the variational parameters. For general graphs this 
step can not be carried on in an exact way, but on locally tree-like graphs it can 
be done with the use of the cavity method on the auxiliary system. For a given 
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instance of the problem, this leads to the following expression of the variational 
energy: 



(id) 



(124) 



where the cavity messages satisfy the usual equations: 

Vivien) = —e b ^ 11 ^^"^fc-rtC *)- ( 125 ) 
k&di\j <? k 

Unfortunately, E v is a complex non-local function of the variational parameters 
Kij and hi, because the cavity messages depend implicitly on far away couplings 
through the cavity equations. 

One way to minimize E v is to introduce the following partition function: 

Z0) = J DiK^Dib^Din^} e -2fe.«*iiM&«M*-*}) /cavity cq , (i 26 ) 

where I C avity cq is a set of delta functions that impose the cavity equations (125) 
and P is a fictitious inverse temperature. In this way, considering the cavity 
messages as independent variables, the computation of Z (more precisely, of 
the average of iV _1 logZ over the disorder, for N — >• oo) can be done through a 
message-passing algorithm whose equations resemble those of the 1RSB classical 
computation. From Z(f3) it is easy to extract the minimum of the variational 
energy by sending the fictitious inverse temperature /3 — > oo. In [231] this 
procedure has been carried out explicitly and reasonable results for the ground 
state energy of the Ising model in transverse field have been obtained; this 
method was then applied to a model of interacting fermions on the Bethe lattice 
in [253], bypassing the usual sign problem. 

5.^.2. Matrix product states 

Within the variational approach, matrix product states (MPS) represent a 
very promising way to study the properties of quantum systems defined on tree- 
like structures. A MPS is a representation of a state of a quantum lattice model 
that is based on a set of tensors defined on the sites of the original model. As we 
will discuss below, this representation is exact on tree-like structures as long as 
the size of the tensors is large enough. A truncation of the tensors size usually 
gives a very good approximation of the state provided the entanglement is not 
too large. 

For one-dimensional systems, MPS are at the basis of many numerical algo- 
rithms, most notably the Density Matrix Rcnormalization Group (DMRG) [254] 
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or the Time Evolving Block Decimation (TEBD) [255]. These methods are 
widely applied in the study of one-dimensional systems where they are known 
to work efficiently while the understanding of their generalization to higher di- 
mensions is still the subject of intense research. 

These methods have been generalized to tree geometries in different works [256, 
257, 258, 259]. The method proposed by [256] is a generalization of the algorithm 
developed by Vidal in [260] in order to study translational invariant systems in 
the thermodynamic limit. For finite trees, related algorithms have been pro- 
posed in [257, 258]. A DMRG algorithm was used in [261] to derive the ground 
state properties of the Hubbard model on the Bethe lattice. 

All these algorithms exploit the MPS representation of the ground state, 
combined with parallel updates that descend from such representation and the 
local properties of the Hamiltonian under study. Once the expression of the 
state in terms of MPS is given one can apply unitary transformations involving 
local operators with update rules that are local [257, 260, 256]. This is a crucial 
property on which the efficiency of the algorithm relies. An important point 
is that during such updates an exact calculation generally brings to increasing 
tensor sizes. However the size of the tensors can be kept fixed exploiting a block 
decimation technique that aims to project on a restricted subspace that carry 
most of the information. 

Beyond unitary operations the same techniques are used to perform the 
imaginary time evolution, which is exploited in the search for the ground state. 
This operation introduces new errors that originate from the normalization of 
the tensors that is spoiled by the non-unitarity of the evolution. Different meth- 
ods have been proposed [256, 257, 260] to account for this effect. 

We refer the reader to the references mentioned above for a more detailed 
discussions of these algorithms. At the end of this section, we will take a slightly 
different perspective by showing how MPS can in principle be used within a 
variational cavity approach. 

MPS for chains. For a chain of N spins with open boundary conditions (a finite 
tree with connectivity 2) a MPS is a state of the form: 

IV') = c CTl ... CTj >i)...kjv> , (127) 

(71,..., fTN 

with 

XI XN 

L 'fTi...(Tjv / j / j lOil Ql 'CK1CK2 CK2 'CK2CK3 ' " - fajV — 1 

a± — 1 ajv — 1 

N Xi 

= II (XNr «*.:)• 

i—l CLi 

where 7M are N matrices defined on the sites of the chain, are A — 1 vectors 
defined on the links of the chain, and in the last equality we set x\^ r = 1, 
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J 1 !"-! _ x J 1 ]"" aTlr ] J"!"" — x J^l 1 " The vpctors A^ arp thp 

Schmidt coefficients that appear in the Schmidt decomposition of the system 
when it is divided into two disjoint parts 1, . . . , i and i+ 1, . . . , N (by "cutting" 
the link between i and i + 1 ) . 

Given a bipartition of the system into two disjoint subparts A(i) = {1, . . . , i} 
and B(i + 1) = + . . . , N}, the Schmidt theorem states that for every vector 
\v) it is possible to find an ortho normal basis for the Hilbcrt space defined over 
A (resp. over B), such that: 

Xi 

K = E A « KU(i)K>s(i+i) , (128) 

a=l 

with = 2 min [ I ' W ~ J ] and Aq > are positive real numbers such that ^* ! =1 (Aq ) 2 
I 

Performing the same decomposition between the sites A{i— 1) = {1, . . . , i— 1} 
and B(i) — {i, . . . , N} one obtains 

Xi-l 

l y ) = E A r 11 MA(i-i) Mb W ■ (129) 

/3=1 

Using the basis defined for the subspace i3(i + 1) one can write: 

Xi 

K)bo) = EE^!^ 1 ^'^ 5 ^ 1 ) ' ( 13 °) 

which defines the matrix 7„ a* used above. In the same way one obtains 

K>A W = E E ^tSi^M^-d • (131) 

CT, ,ti = l 

The orthonormality of the basis |i> a )yt(i) and |f a )_B(i) imposes normalization 
conditions on the A's and 7's. They must indeed satisfy, Vi = 1, . . . , N: 

E(aW) 2 = i, 

Q = l 

Xi 

B®{v a '\v a ) B (i) = E E a J' (tS)* ^J 1 = 5 a , a > , (132) 

a-j ,3=1 
Xi-l 

/ I \ \ " \ " [i]ffi>[i-l]/ [i]<Ti\* \ [i-X] c 

<Ti = 1 

The average values of local observables can be easily computed. Let us consider 
a one-body operator 

° W = E°?,J^|. (133) 
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Then 



WOWW = E alI ^0^ 4 E^E^AP(7a^rAMA^ 1] (7a?)AL fl , 
while for a product of two one-body operators 

w w= 2 o^oSS^I^X^^I ( 134 ) 

(Ti .(J^ .a'- 

the expectation value is: 

(135) 

The expression (127) is exact if the dimensions Xi 01 the A's and T's are 
large enough, however the same definition can be used in a variational way, 
for fixed (small) sizes, and still provides a good representation of the ground 
state. The entanglement is usually used as a measure of the accuracy of such 
representation. The Schmidt coefficients in fact are directly related to the so- 
called entanglement entropy through the formula: 

Xi 

S A = -Tr [p A log Pa] = -J2 ( A « ) 2 lo g[( A a ) 2 ] = S * . PA <* Tr B . 

a=l 

(136) 

In the limit in which Xi = lj i- e - if the two parts of system are separable 
\v) = \v)A(i)\v)B(i+i), then S A = 0. 

MPS for trees. Trees have no loops, thus, removing an edge divides the system 
into disjoint parts. The Schmidt decomposition can be applied and it allows to 
naturally define MPS also in this context. We refer the reader to [256, 257, 258] 
for more details and state in the following the expressions derived in these works. 

In the case of trees the expression (127) is generalized using a vector A« for 
each edge and tensors "fa\ a !...a c with c lower indices (where c is the connectivity of 
the graph) plus one spin index as in the one-dimensional case. The normalization 
conditions generalize for all the c indices of the tensors. In order to derive the 
vector Xa one has to perform the Schmidt decomposition on the corresponding 
edge (ij). This divides the system into two disjoint subtrees 7l->j and Tj^i, 
where each subset contains the connected component made of sites connected 
to i and j respectively: 

x<y> 

\v)=J2x^\v a ) Ti ^\v a ) T ^ i ■ (137) 
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The tensors jau'..,a c are obtained performing the Schmidt decomposition for 
the c bonds surrounding the site i and then expressing one of the orthonormal 
basis that derive in terms of basis obtained with the other c — 1 decomposition 
and the spin similarly to the one-dimensional case: 

k,,)^ = E E 7{?U 6<N n [ttw^w- 

<Tj {a (lfc) } : fceai\i fce9i\j 

(138) 

where each cxm^ is summed from 1 to XUk) > with normalization conditions 
completely analogous to Eq. (132). In this way one arrives at the following form 
for a general vector, in terms of matrix product states: 

N 

m=e e (n^u 6 J(n«,)i^>---i^>- (139) 

{<r»} {<*<«>} <«fe) 

Expectation values of one-body operators are given by 

mo%> = E °EU E <:, }fce8 J<L >} _r n [«>] 2 ■ 

In order to compute the expectation value of a two-body operator O^O^l acting 
on two sites i and j we denote with S the set of sites on the unique path joining i 
and j, with V the edges that are adjacent to at least one vertex in SU{i,j}, and 
with 1Z the set of the edges that are adjacent to exactly one vertex in <SU 
More explicitly, V\TZ are the edges of the path between i and j. Then one has 



x v ^ h m r-y b1<7j r-v^f* i* 

'{a(ik>}fceai L '{a' (ifc) }fceai J '{ Q (jfc> Ifceaj 1 i a '(jk) Heaj 1 

{a (lk) },{a' <lk) }:{lk)eV 

V TT^W^ lW'i 1* TT >< ife ) x^ fc > TT A 

ll 7 {a<ife)}feeei L7 {«' <!fe) }feeei J 11 11 ° Q ("=> . Q < ifc) ' 

The disordered Ising model in transverse field. To show why MPS are particu- 
larly useful, we can consider again the simplest example of an Ising model on 
a regular graph of connectivity c, with Hamiltonian H = — ^2(ij) Jij^i^j ~ 
Y^i^i&fi an< i a variational MPS \ijj v ) with tensors of fixed size \. Then the 
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variational energy is 
E v = (MH\^) = -Y.^ E 7fc[7& a ]*II(^ > ) a 

-E^E-' E f n (a^) 2 ) 7£ }5 (7{t>, 



XA 5 A <5' 7,5 {/ 5 fc} l7 5{/ 3 fc }j ( 11 l*0 k ) 



(140) 



with normalization conditions (on each directed link): 

£(W = i, 

/ \ (141) 

E E II ^) 2 7[t } ,(7{t^r=^. 

CT {<*k}kedi\j \k£di\j ) 

If we consider first a model without disorder, i.e. with = J on all edges 
of an infinite tree, and r, = V on all vertices, then all sites are equivalent and 
we can assume that the tensors and vectors do not depend on the site and link 
indices. The variational energy is then a function of a finite set of variational 
parameters and one can devise several strategies to minimize it, either based on 
numerical minimization routines, or on simulated annealing. Alternatively, one 
can use the strategy of [256, 257, 258] by applying the imaginary time evolution 
to the variational state. This procedure gives the variational energy in the 
thermodynamic limit. 

Cavity optimization of MPS. In the case of a disordered model with random- 
ness in the couplings and/or local transverse fields one should keep all vectors 
A and tensors 7 as variational parameters in the expression (140). As a direct 
comparison of Eq. (140) and Eq. (124) shows very explicitly, there is a crucial 
advantage of MPS with respect to e.g. Jastrow wavef unctions. For MPS, the 
variational energy can be written as an explicit function of the variational pa- 
rameters, which is a sum of local terms involving the tensors on a site i and its 
neighbors and the vectors on the links among these sites. On the contrary, for 
the Jastrow wavefunction a cavity computation is needed to write E v , which is 
therefore a very implicit expression of the variational parameters. 

For a generic MPS describing a disordered system, the variational energy can 
then be interpreted as a "classical Hamiltonian" for a system whose classical 
variables are the tensors and vectors of the MPS. The partition function of such 
a model would read schematically as 

20) = f d{7 [ll }4A <lJ> }e-^ [ ^ W} ' {A<IJ)}I /„or m aiization[{7 W },{A fo> }] (142) 
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with E v given in Eq. (140), /normalization a set of delta functions enforcing the 
normalizations in Eq. (141), and ft a fictitious inverse temperature to be sent 
to oo at the end. 

The latter is much simpler than Eq. (126) where a functional delta over the 
cavity messages appears; this is not required here because the variational energy 
for MPS is an explicit function of the tensors. The partition function (142) can 
then in principle be computed via the standard classical cavity method 6 . The 
price to pay is of course that the basic classical variables are tensors of size x> an d 
the cavity messages are therefore distributions over the space of such matrices. 
The limit j3 — > oo can be performed explicitly in this case following [153], leading 
to the optimized variational energy. Although this strategy was not turned into 
a concrete calculation for the moment, we believe that it is a very promising way 
to compute the zero temperature properties of quantum random optimization 
problems, and therefore complementary to the finite temperature PIQC. 

5.5. Exact diagonalization and numerical integration of the Schrodinger equa- 
tion 

We now present briefly exact numerical techniques to study the statics and 
the dynamics of finite quantum systems. If the size of the Hilbert space is small 
enough, thermodynamic and spectral properties of H can be obtained from 
exact diagonalization techniques. Such techniques are a whole field of research 
in themselves so we just give a brief overview here. If one is interested in finding 
all the eigenvectors and eigenvalues of a given matrix, the most commonly used 
techniques are the Jacobi and the Gauss-Seidel methods. If one is interested only 
in the low energy part of the spectrum, it is possible to use Lanczos type methods 
to obtain more quickly the lowest lying eigenvectors. Moreover, these methods 
require only to be able to compute the multiplication of a vector of the Hilbert 
space by H . In the case of sparse matrices as the ones relevant for optimization 
problems, this can be done without keeping the whole matrix H in memory; 
therefore the limitation of this method comes from the size d N of a vector of 
the Hilbert space (where d is the dimension of a single qudit), rather than from 
the size d 2N of the Hamiltonian operator. For quantum 1/2 spins where d = 2, 
this allows for computations up to TV = 25 on standard computers. The results 
reported in Sec. 6 were obtained by using the Arpack package [262]. Note that 
in some cases, the presence of exact symmetries (operators that commute with 
H) allows to reduce the size of the Hilbert space and thus increase the sizes of 
the systems that can be exactly diagonalizcd. 

For the analysis of the quantum adiabatic algorithm it is particularly inter- 
esting to simulate exactly the real time evolution of a quantum system following 



6 The presence of interactions involving a variable and all of its neighbors requires using a 
"trick" consisting in creating a copy of each variable on its neighboring sites: sec [198] for a 
more detailed discussion of this point. 
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the time-dependent Schrodinger equation denned in Eq. (12): 

lA Ws))= £( sMs)) . (143) 

This is a linear differential equation, which can thus be solved using standard 
numerical integration techniques, such as Runge-Kutta or Adams-Bashforth 
methods. However, it is better in practice to make use of the Hermitian nature 
of the generator of the dynamics [263]. In fact, (143) can be rewritten as: 

= T (e- lT 'o |^(o)) = U(Q, sM(Q)) (144) 

where T denotes the time-ordering operator. Because H(s) is Hermitian, U(0, s) 
is unitary; one is then interested in finding unitary approximations to U(0, 1). 
The simplest way to do it is to write: 



U 



71 11 

(0, 1) = T (e- iT /o 5 ^ ds ) = JJ T (e- lT ^ +1 R ^ ds ) = \{U{ Sl , s i+1 ) 



i=l i=l 

(145) 

We are interested in the particular case of a linear dependency of H(s) on s: 
H(s) = (1 — s)Hi + sHf. The approximation 

iT/; +Aa s'Hids'\ ( p -iT /;+ As (l-s')i?ids' 

(146) 
= U As {s) 

gives rise to an error in operator norm \\A\\ = sup x \\x=i\\ ll^^l bounded 
by [264, 265]: 

\\U(s,s + As)-U As (s)\\ <\\[H u Hi]\\T^L + 0(As 3 )=0(NTAs 2 ) . (147) 

Indeed in all the cases of interest here the commutator of the initial and final 
Hamiltonian has a norm ^of order N. We define the approximate evolution 
operator U(0, Sj) = Y[]=o U As (sj). The triangle inequality 

\\U(0,s i+1 )-U(0,s i+1 )\\ = 

= \\U(0,8i)(U(8i,8i +1 ) - U As ( Si )) + {U{0, Si ) - U(0,Si))U As ( Si )\\ (148) 
< \\U( Si , Si + As) - U As (si)\\ + \\U(0,Si) - U(0, Si )\\ 

leads by recurrence to 

\\U(0, 1) -U(0, 1)|| < 0{nNTAs 2 ) = 0{NT/n) . (149) 
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One can thus replace the exact evolution operator U(0, 1) by its approximation 
C7(0, 1) with a precision of order e in the evaluation of intensive observablcs if 
the number of discretization steps n is of order NT/e. Note that this bound 
does not involve the spectral gap of the system, which is thus not directly the 
bottleneck for the simulation of the quantum evolution. 

Let us evaluate the total complexity of the procedure. Doing the computa- 
tion in Hf eigenbasis, the multiplication by an operator e aHf is trivial and can 
be realized in a time proportional to d N . where d is the dimension of the Hilbert 
space of a single qudit. On the other hand, taking for H\ a sum of identical op- 
erators acting on single qudits, Hi = Y^,j=i fy/j wg can write e aHi = Yif=i s ahj ; 
the exponential one has to compute is the same for any site and its action on 
a vector of the Hilbert space can be computed with less than Nd 2 operations. 
Therefore we finally see that the action of U& s {s) on a vector can be computed 
within 0(d N ) operations; leading finally to a number of operations bounded by 
NTd N /e to obtain a precision e on the final result of the evolution. Practically, 
the resources limitations of this method come both from the size d N of the vec- 
tor one has to keep in memory, and from the large time T one is interested in 
for quantum adiabatic computations. 

As a side-remark let us mention that the matrix product states approximate 
paramctrization of quantum vectors, discussed in Sec. 5.4.2, can also be used to 
study the real-time (Schrodingcr) dynamics of quantum systems, see [266, 267, 
268] for details. 

5.6. Quantum Monte Carlo 

In this section we discuss how Quantum Monte Carlo (QMC) simulation 
algorithms can be used to extract relevant information on random optimization 
problems, in particular their energy gap. 

Path Integral Quantum Monte Carlo (PIMC) simulations have a very long 
history and were initially performed to study continuum systems with particular 
focus on Helium 4 [269]. Rapidly, several implementations were developed to 
study lattice systems, made of spin or bosonic degrees of freedom. Early imple- 
mentations were based on a Suzuki- Trotter path integral in discrete imaginary 
time [270, 271, 272], but rapidly it was realized that the continuum imaginary 
time limit could be taken explicitly [273, 274, 275, 224, 276]. A similar approach 
is based on an exact sampling of the pcrturbative expansion and goes under the 
name of Stochastic Series Expansion (SSE) [277, 278]. 

PIMC is a Monte Carlo method that produces configurations of imaginary 
time trajectories sampled from the measure fiQ we introduced in Eq. (93). In 
the case of spins 1/2 in a transverse field, i.e. for H = J2 a E(q)\g) (<x| — V J^. erf , 
it reads more explicitly: 

i N 

= ^(0),^)ll rkil e-tf , (150) 
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with 

r N 

Z = Tre-? H = / JjD < r i rl«"l e-/o"«(2(*))d* , (151) 

z(0)=z(/3) i=1 

where |o^| is the number of flips in the trajectory <7j(i) of the i'th spin. Empiri- 
cal averages over the trajectory configurations allows to compute the thermody- 
namic (both thermal and quantum) averages (•) = Tr [•e~@ H ]/Z. The various 
versions of PIMC differ in the allowed moves between configurations a_ , that are 
usually required to fulfill the detailed balance condition with respect to the mea- 
sure /xq, to ensure the stationarity of the latter. The PIMC results reported in 
Sec. 6 have been obtained by using the heat-bath algorithm introduced in [224], 
in which the PIMC updates consist in drawing a new trajectory for a randomly 
chosen spin, according to its conditional probability induced by its neighbors. 
This is possible thanks to the path generation procedure explained in Sec. 5.2.3; 
a rigorous proof of fast convergence for this algorithm can be found in [246] for 
the Ising ferromagnetic model in transverse field on an infinite tree. 

As explained in Sec. 2.3.2 the efficiency of the quantum adiabatic algorithm 
is controlled by the gap between the two lowest lying eigenstates of the interpo- 
lating Hamiltonian; in the following wc discuss two different strategics to extract 
such a gap from QMC simulations. 

5.6.1. Extracting the gap from correlation functions 

The first strategy [279, 197, 278, 280] is based on the computation of imag- 
inary time correlations. Consider for example the spin-spin correlation 

(a?(r)5?(0)) = iTr [e-C-^Sfe-^o?] , (152 ) 

that as any other observable can be easily computed via PIMC. Let us denote 
by Eo < Ei < . . . the distinct eigenvalues of H, with associated eigenvectors 
\n). In the limit (3 — > oo with r fixed, inserting the representation of the identity 
/ = \n)(n\, one obtains the spectral representation of this function as 

(3f (r)3f (0)> = £ \{0\d!\n)\ 2 e-^- E °) . (153) 

n 

Then, if the limit r — > oo is taken (after /3 oo), we have 

(5f(r)Sf(0))-(5f> 2 ~e- rA , (154) 

where we denoted A = E\ — Eq the energy gap between the two lowest levels 
(this formula is easily generalized if the ground state and first excited level are 
degenerate). Even though PIMC simulations can only be performed at finite 
f3, in the regime 1/A <C r <C j3 the plot of the logarithm of (of (r)af (0)) c 
versus r is a straight line that can be fitted to extract A. We refer the reader 
to [279, 197, 278, 280] for details and concrete examples of such computations, 
in particular to [280] where an optimal choice of observables in the computed 
correlation function is discussed. 
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5.6.2. Extracting the gap from the specific heat 

Another possible strategy to compute the energy gap A is based on the 
evaluation of the specific heat: 



C= ±(H)=^((H 2 )-(Hf) 



(155) 



If again E < E\ < 
degeneracies go, g%, . . . 
heat behaves as 



. . are the distinct eigenvalues of H, with associated 
one sees that in the low temperature limit the specific 



C 



2 91 PA 



(f3A)--e 

go 



(156) 



The value of A can thus be obtained from the behavior of the specific heat at 
low temperatures. Moreover C can be computed from a QMC simulation. Let 
us introduce a fictitious parameter x and define Z(x) = Tr e~@ xH , in such a way 
that 



(H) 



1 Z'(x) 



13 Z{x) 



(H 



x=l 



1 Z"{x) 
~fP Z(x) 



As we explained above, with the path integral representation we have 

N 



Z(x) = 



(157) 



(158) 



On this form it is very easy to take the derivatives with respect to x, which 
leads to 



(H) = 



(H 2 ) 



Da^ia) | ~ jT" E(a(t))dt ~ r E ^ } 
1 r0 N l l' 



(159) 



1 N 

' i=l 



These two quantities, and in consequence the specific heat, can thus be directly 
determined from configurations of paths generated in a QMC simulation. The 
advantage of this procedure with respect to the previous one is that in prin- 
ciple the specific heat should be easier to compute than the time-dependent 
correlation functions, and that there is no need to identify the correct regime 
1/A <C r < ,8 in r. The disadvantage is however that one has to perform sim- 
ulations at several temperatures (3 > 1/A to extract the slope of log(C) plotted 
as a function of /?. 



5.6.3. Imaginary time annealing 

Let us finally mention two numerical approaches that, although very differ- 
ent, can be both termed "imaginary time annealings" . 

A first "imaginary time annealing" consists in solving the Schrodinger equa- 
tion in imaginary time (which can be done either by exact diagonalization or 
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by QMC), i.e. to study the evolution of a vector of the Hilbert space \tp{s)) 
according to 

-J^M')) = H(s)\tP(s)) , H(s) = (1 - a)Hi + sHf , (160) 

to be compared with Eq. (12). This evolution is not unitary, hence would not 
be realizable with a quantum computer, but can be implemented numerically. 
The reader will find in [21] and references therein a detailed comparison of the 
real and imaginary time version of the quantum annealing. 

A distinct procedure, that we shall use in Sec. 6, is more precisely an an- 
nealing of the Path Integral Monte Carlo procedure. In PIMC, configurations 
of imaginary time paths are generated with the measure (1q of Eq. (93), by 
performing many Monte Carlo updates on the configuration to approach this 
stationary measure. In other words, as in any Monte Carlo procedure, one starts 
with a given initial configuration, and repeatedly apply to it a given operation 
to produce new configurations. One therefore introduces a fictitious time tyic 
that describes the number of such operations that were done since the beginning 
of the procedure. In a PIMC annealing, the parameters of the measure (inverse 
temperature (3 and transverse field T) slowly evolve during the PIMC simula- 
tion, and become "Monte Carlo time" -dependent parameters /3(£mc)j r(^Mc)- 
Note that classical simulated annealing is a particular case of this procedure 
where T = at all Monte Carlo times; on the other hand one can set (3 to a 
very large fixed value (very small temperature) and let T evolve. In the limit 
where /3 is infinite, and the rate of variation of T vanishes, this coincides with 
an adiabatic Schrodinger evolution: at all Monte Carlo times the configuration 
of paths is drawn from the measure (iq which encodes the instantaneous ground 
state of the original quantum Hamiltonian. Therefore, PIMC annealing allows 
for an interesting interpolation between classical simulated annealing and zero 
temperature quantum annealing. Note however that the condition of adiabatic- 
ity for the PIMC annealing has a priori nothing to do with the one of the original 
Schrodinger evolution (see however [281] and references therein for a discussion 
of this point). The relevant gap is in the latter case the one of the quantum 
Hamiltonian H, while in the former case it is the one of the Fokker-Planck 
generator of the PIMC dynamics on the space of path configurations; we refer 
to [24, 246] for further analysis of this PIMC annealing. As a final important 
remark, note that the clustering transition (or "dynamic transition"), which is 
signaled by the appearance of a non-trivial solution of the 1RSB cavity equa- 
tions at m = 1 (both in the classical and quantum cavity method), is directly 
related to a glassy lack of equilibration of the PIMC (if the thermodynamic 
limit is taken before the limit of infinitely slow annealing) [179]. The decorre- 
lation time of PIMC dynamics in tr^c becomes infinite at this transition. This 
provides a very useful way to detect this transition using PIMC, that we will 
illustrate in Sec. 6. 
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6. Results on specific random optimization problems 

The aim of this section is to apply the methods outlined in Sec. 5 to the 
study of random instances of real optimization problems, such as those defined 
in Sec. 2.1.1. At variance with the "toy" models investigated in Sec. 4, the 
problems that we will discuss in this section are standard problems in computer 
science. Still we expect that their phenomenology is close to the one of the 
toy models. An important lesson that we learned from Sec. 3 and Sec. 4 is 
that there is a wide variety of behaviors in the different classical optimization 
problems, that only become more complicated when quantum fluctuations are 
added. In consequence we shall not try here to propose a complete classification, 
but present results on a few specific models (or "case studies"), that illustrate 
the main phenomena that were discussed in Sec. 4 and lead to exponentially 
small gaps, namely first order transitions and level crossings. Part of these 
results were already published, part are original. The analysis is based on the 
methods that have been described in Sec. 5, we concentrate here on the physical 
results. 

6.1. Early results 

To put the discussion in a historical perspective, it is worth to mention that 
most of the current interest in the Quantum Adiabatic Algorithm (QAA) was 
triggered by a series of early numerical works that found evidence for a polyno- 
mial scaling of the minimum gap in the l-in-3 SAT (or Exact Cover) problem 
(see Sec. 2.1.1) with a transverse field, suggesting an exponential speedup with 
respect to classical algorithms [20, 279]. These studies considered a particular 
ensemble of random Exact Cover instances, constructed to have a Unique Sat- 
isfying Assignment (USA); see [282] for a detailed description of the procedure. 
We call this ensemble EC-USA in the following. This choice was made because 
in this case the minimal gap between the ground state and the first excited state 
can be unambiguously defined. 

The original work by Farhi et al. [20] was based on exact diagonalization 
of very small (N < 20) EC-USA instances. A polynomial scaling of the gap 
for those instances was detected. This was initially confirmed by a Quantum 
Monte Carlo (QMC) study of EC-USA instances with N < 128 [279]. However, 
it was shown later by the same authors [197], by using a more refined analysis 
and larger (N < 256) sizes, that an increasing (with N) number of instances 
display a first order transition and an exponentially small gap. 

These numerical studies were extremely difficult not only because of the 
numerical cost of the exact diagonalization step (Sec. 5.5), but also because 
EC-USA instances have an exponentially small probability in the fully random 
ensemble of Exact Cover [282]. Therefore, already constructing EC-USA in- 
stances is an exponentially hard task (Sec. 3.9) and constitutes one of the main 
limitations to access large sizes [20, 279, 197]. Furthermore, strong finite size 
effects were detected on EC-USA instances [197]. 

Another important problem, that has been discussed in more details in 
Sec. 3.9, is the following. Suppose that we take a fully random ensemble of in- 
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stances of a given problem, that are exponentially hard for known classical algo- 
rithms with probability one when N — > oo. If one selects only the USA instances 
of this ensemble, and if the latter have exponentially small probability, then one 
is conditioning on extremely rare instances and it is not obvious anymore that 
these instances remain exponentially hard for classical algorithms [196]. There- 
fore, despite some numerical evidence for an exponential scaling of the running 
time of classical algorithms for EC-USA of N < 256 [279, 282] has been re- 
ported, the classical computational complexity of EC-USA at much larger N 
remains an open problem, even if rather academic as the instances of this prob- 
lem essentially do not exist in the thermodynamic limit. 

These early studies highlighted the importance of being able to construct 
USA instances of much larger sizes (Sec. 3.9), and of being able to investigate 
such instances in the thermodynamic limit via the quantum cavity methods 
(Sec. 5). Both these goals can be achieved by using the so-called "locked mod- 
els" [191]. We consider one of such models in the next section. 

6.2. Locked models: XORSAT on a regular graph 

In this section, we examine the XORSAT problem on a random regular 
graph, a typical representative of the class of locked models [191] that were 
discussed in Sec. 3.9. The main property of these models is that clusters of 
ground state configurations do not have internal entropy: they are isolated 
points. Therefore we do not expect level crossings induced by the energy-entropy 
competition discussed in Sec. 4, which simplifies a lot the analysis of the mod- 
els. Moreover, these are the simplest models to study with the cavity method, 
allowing us to illustrate the usefulness of the method in the simplest non-trivial 
setting. A crucial properties of these models is that the parameters can be tuned 
in such a way that USA instances have a finite probability for JV-fooin the 
random ensemble, see Sec. 3.9. 

Preliminary results we obtained on the quantum XORSAT model were re- 
ported in [48] , where we showed the existence of a first order transition associ- 
ated to an exponentially small gap in N. We present these results in much more 
detail in the rest of this section, together with some previously unpublished re- 
sults. Other locked models have been studied in [278, 283] and showed a similar 
behavior. 

6.2.1. Definition of the model and its classical properties 

We focus on the fc-XORSAT problem, defined on a random c-regular graph, 
which has been studied in the classical case in [284, 285], to which we add a 
quantum transverse field. In quantum spin language, the model is defined by 
the following Hamiltonian (where we omit a factor of 2 with respect to the 
definition of Sec. 2.1.1): 

M N 
H = H P + TH Q = ^2(1- J SJ....SJ.)-r^3f . (161) 

a—1 i— 1 



in 



Here, J a = ±1 with equal probability. The k spins i?,^, ••• ^1 involved in 
clauses a = 1, ■ ■ ■ , M = Nc/k are chosen uniformly at random among all possi- 
ble choices such that each spin enters exactly in c clauses. This defines a regular 
random graph structure where variables have connectivity c and interactions 
have connectivity k. 

As usual, in the classical limit T = 0, a given instance of the problem (defined 
by the choice of the random graph and of the couplings J a ) is called satisfiable 
(SAT) if there is a ground state of zero energy, UNSAT otherwise. It is easy to 
see that the annealed entropy (i.e. the logarithm of the average number of so- 
lutions) density is log(2) (1 — c/k) when T = 0. It has been shown in [284, 285] 
that when the annealed entropy is positive (c < k) the model is SAT with a 
probability going to 1 as N — > oo and the typical number of solutions is exponen- 
tial in N, concentrated around its mean predicted by the annealed entropy. On 
the contrary if the annealed entropy is negative (for c > k) satisfiable instances 
are exponentially rare, and typically the model is UNSAT. In the marginal case 
c = k the model is SAT with finite probability, and when it is SAT the number 
of solutions is typically finite. 

We are particularly interested in USA instances of Hp. Based on the above 
discussion, it is clear that these instances are exponentially rare if c ^ k, and 
it is natural to expect that they have a finite probability for c = k. We have 
indeed found numerically that for c = k = 3, in the limit N— >oo, the fraction 
of SAT and USA instances are / SA t = 0.609 ± 0.003 and / USA = 0.2850 ± 
0.0022, as determined by using either Gaussian elimination, or a Davis-Putnam- 
Logemann-Loveland-like algorithm [286], to count the number of solutions of 
40000 instances of different sizes and extrapolating the result to N — > oo [48]. 
Similar values are obtained for c = k = 5 and c = k = 7 (for even values of k 
solutions always come in pairs due to the spin flip symmetry). It is worth to 
mention that in the limit k = c — >• oo (taken after N —> oo), USA instances 
of the model should approach a particular Quantum Random Energy Model 
(QREM) where the distribution of the classical energies is a binomial. This 
model was analyzed in [41], and the existence of a first order transition was 
shown rigorously, supporting the results obtained with the cavity method at 
finite k and c. See [282] for numerical studies of /us A in other locked models. 

For this model, Path Integral Quantum Cavity (PIQC) computations have 
been performed following the method described in Sec. 5.2. The model is fac- 
torized (Sec. 5.1.5) thanks to the regular structure of the hypergraph, hence 
the replica symmetric (RS) computation requires a single population of A/traj 
imaginary time trajectories (we used A/traj = 10 5 ). For 1-step replica symmetry 
breaking (1RSB) computations, we used A/int =4000 populations of A/traj =4000 
trajectories. In both cases, each data point has been obtained as an average 
over 100 cavity iterations. We checked that finite population size effects are 
negligible with this choice of parameters. In addition, we will report the behav- 
ior of the spectral gap, determined by means of Exact Diagonalization (ED), 
and some Path Integral Quantum Monte Carlo (PIMC) results. 
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Figure 18: Phase diagram of the fc-XORSAT model Eq. (161) on a e-regular random graph. 
The top panel represents the (T, T) phase diagram, displaying four possible phases: CP, dCP, 
QP, SG (see text). Open symbols are RS results: first order transition line Tf D (r) (open 
circles) separating the CP and QP, with the corresponding spinodals (open diamonds). Full 
symbols are 1RSB results: clustering transition 7d(r) (full squares) separating the CP and 
dCP, condensation transition T C (T) (full triangles) separating dCP and SG, T{ a (T) (full circles) 
separating the SG and QP. The middle panel reports the transverse magnetization m x as a 
function of T, and the bottom panel reports the free energy density from PIQC or the energy 
density from PIMC as a function of T, all at fixed temperature T = 0.05. In these panels, full 
lines are PIQC results (RS or 1RSB) while symbols are PIMC results. 

(Left panel) k = 4 and c = 3. There is no SG phase. The CP phase becomes a dCP at low 
enough temperature, while a first order transition separates the CP (or dCP) and QP phases. 
m x jumps at the first order transition. 

(Center panel) k = c = 3. Also here there is no SG phase. PIMC data are reported, for 
a sample with N = 2049: red diamonds are obtained starting from the QP (r = 2) and 
decreasing T, while black squares arc obtained starting from a classical ground state (found 
using Gaussian elimination) and increasing T. 

(Right panel) k = 3 and c = 4. Here a SG phase is present, and Tf (r) separates the SG and 
QP phases. PIMC data are reported for N = 120 and averaged over 20 samples (full symbols) 
and extrapolated in l/N to the N—too limit (open symbols). Black curve, starting from the 
classical ground state found using an exact MAXSAT solver [287]. Red curve, starting from 
the QP. 
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Figure 19: Complexity of the 4-XORSAT model at c = 3. (Main panel) Equilibrium complex- 
ity as a function of T at fixed T. The complexity remains finite up to the dynamic transition 
line Td(T), where it jumps abruptly to zero. (Inset) Equilibrium complexity of the classical 
model (T = 0) as a function of temperature. The dashed line marks the T = value, log(2)/4. 



6.2.2. Exponentially degenerate ground state: c < k 

As a representative of the case c < k, we take here k = 4 > c = 3. The clas- 
sical ground state is exponentially degenerate with entropy log(2) (1 — c/k) = 
log(2)/4. It can be shown via the cavity and replica methods [284, 285, 191] or 
using rigorous methods [167, 168] that the ground states are arranged in isolated 
clusters. Therefore, the internal entropy of each cluster is zero, the complexity 
of clusters is E = log(2)/4, and typically the clusters (solutions) have Hamming 
distance of order N, therefore they are very far away in configuration space. 
The classical equilibrium complexity as a function of temperature is plotted in 
the inset of Fig. 19. The model is SAT with probability one and the typical 
number of ground states is exp(AE) = 2 N ' 4 . 

The phase diagram of the model, as obtained from PIQC, is reported in 
Fig. 18. The RS computation predicts, at low enough temperature T < 0.3, 
a first order transition between two different paramagnetic (m z = (erf) = 0) 
phases: the Classical Paramagnet (CP) characterized by a small value of trans- 
verse magnetization m x = (erf), and the Quantum Paramagnet (QP) that has 
a larger value of m x . The first order transition is signaled by a jump of m x 
and a crossing of the free energies of the two phases, that can be clearly seen in 
Fig. 18. As in any first order transition for a mean field model, the two phases 
can be continued in the region where they are metastable until a well defined 
spinodal point. The transition line and the corresponding spinodals are shown 
in the {T,T) phase diagram in Fig. 18; the transition is found at a slightly 
temperature-dependent Tf (T) « 3/4 = c/k. We also report in Fig. 18 the cav- 
ity method predictions for m x and the free energy density / = —T\og(Z)/N at 
very low temperature T = 0.05. 

Next, we discuss the outcome of the 1RSB computation. As discussed in 
Sec. 5, the key quantity that is computed in this approach is the equilibrium 
complexity E eq (r,T) (at m = 1), which is reported in Fig. 19 as a function 
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of T and T. Below the classical dynamical transition T^iY = 0) ~ 0.41, the 
complexity is positive in the classical case (see the inset of Fig. 19). Increasing 
r, we found that E eq (r,T) remains independent of T, until the dynamic tran- 
sition line Td(T) is met, then the 1RSB solution disappears discontinuously, as 
revealed by the abrupt jump of the complexity curves in Fig. 19. The fact that 
E cq is independent of Y is not surprising, based on the discussion of the random 
subcubes model (Sec. 3.6). If the clusters do not have any internal entropy, 
and if their relative Hamming distance is of order N, different solutions are not 
mixed at any finite order of perturbation theory in Y. Therefore, each classical 
ground state is continuously transformed in a quantum eigenstate. Moreover, 
the local environment around each ground state is the same for N — > oo, hence 
at any finite order of perturbation theory the quantum energy is the same for all 
ground states, and the degeneracy is not lifted. This can be shown as follows: 
suppose that r = {r^} is an assignment of the classical spins corresponding to 
a ground state. Then J Tji • • • = 1, Va. We can apply to the Hamiltonian 
(161) the unitary transformation U T = {af — > T^crf } and we obtain 

M N 

£ = ..s%)-r£sf . (162) 

a=l i=l 

Therefore, for each classical ground state r there is a symmetry that allows to 
map the Hamiltonian into a ferromagnetic one and map r to the ferromagnetic 
state. This shows that perturbation theory around each classical ground state 
gives identical results. The number of ground states remains constant and equal 
to its classical value, 2 N ' i , so that the complexity is constant as a function of 
r. Moreover, the fact that the complexity remains approximately constant at 
all temperatures suggests that there are no level crossings between states of 
different intensive energy, as in the QREM [38] . 

The main result is then that the equilibrium complexity at to = 1 is positive 
or zero everywhere. The implications of this result are twofold: first of all, it 
confirms that the RS computation of the thermodynamic observables is in this 
case correct in the whole phase diagram (r, T) (remember that as discussed 
in Sec. 5, a true 1RSB phase is signaled by a negative complexity at to = 1). 
Therefore, the only thermodynamic singularity is on the first order RS transition 
line. Secondly, the complexity is strictly positive for low enough values of T and 
r, implying that the CP phase is actually a dynamical CP (dCP, the meaning 
and motivation for this name have been discussed in Sec. 5) where an exponential 
number of states coexist. The point where the equilibrium complexity becomes 
positive is the clustering temperature T&(Y), and is reported in Fig. 18. Wc 
recall that equilibrium thermodynamic properties are unaffected as one crosses 
the transition between CP and dCP (Sec. 5.1.4). 

In Fig. 20 we show ED results for this case. The lowest part of the spectrum 
of a typical instance with N = 16 is plotted as a function of Y. The instance 
we show has a ground state degeneracy N = 2 N / A = 16 at Y = 0, which is 
the most probable value. Increasing Y, we see that the lowest 16 levels remain 
extremely close in energy (the difference is expected to be exponentially small), 
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up to a value of T ps 0.75, the location of the hrst order transition in the 
thermodynamic limit. At this value of T we observe that the 17th state (the 
first classical excited state) goes down in energy and approaches the bunch of 
ground states. The figure suggests the presence of an avoided crossing between 
this state and the set of ground states. These data confirm the cavity prediction: 
the ground state remains exponentially degenerate at any finite T < Tf Q , while 
at Tf a first order transition happens, caused by a level crossing between the 
degenerate pure states of the dCP and the QP. 

The determination of the gap is complicated by the fact that for a given 
instance the ground state has degeneracy Af, the average of N~ 1 \ogAf over 
instances being equal to the zero temperature complexity log(2)/4. Therefore, 
the interesting gap to determine the performances of QAA is the minimal gap 
(over T) between the lowest energy state and the (Af + l)-th excited state, that 
we call A m i n : transitions to any lower energy state are not dangerous because 
these states will continuously transform into one of the classically degenerate 
ground states. A more detailed discussion of the QAA dynamics in presence of 
almost degenerate levels might be in order here but, as discussed in Sec. 2.3.2, 
no precise results are available at present. Because Af increases exponentially 
fast in iV (it concentrates quickly around the value 2 N / 4 ), one has to compute an 
exponentially large number of levels, which slows down considerably the exact 
diagonalization code 7 and in practice limits us to N < 20. In Fig. 21 we report 
data for the minimum gap as defined above. Despite the strong size limitations, 
the scaling of the gap appears to be exponential in N, as expected at a first order 
transition. We observed that fluctuations in Af induce large fluctuations of the 
gap: indeed, restricting the average to instances having exactly 2 N / 4 classical 
ground states reduces a lot the fluctuations, and the curve is much closer to an 
exponential, at the same time the difference at large ./V being extremely small 
(we don't show the corresponding data). We will discuss further the behavior 
of the gap at the transition in the simpler case c = k. which we analyze next. 

6.2.3. Finitely degenerate ground state: c—k 

We now turn to the case c = fc, where the complexity at T — vanishes 
and the number of ground states is finite with finite probability. We choose 
as the simplest example c = k = 3. The phase diagram, reported in Fig. 18, 
is qualitatively identical to the one we obtained for c < fc, the only difference 
being that the equilibrium complexity now vanishes for T = and any T > 0, so 
the number of ground states is finite for any T. Another quantitative difference 



7 This problem could be removed by making use of the symmetries U T discussed above. It is 
reasonable to assume that both the ground state and the (A/"+l)-th excited state belong to the 
subspacc of the Hilbert space that is completely symmetric under all the J\f such symmetries. 
Then one could restrict the Hamiltonian to this subspace and compute the gap between the 
two lowest levels. The resulting reduced matrix will not be sparse, so this strategy is not 
efficient for exact diagonalization and it was not used here. However, a similar strategy turns 
out to be very useful if one wants to compute the relevant gap via QMC, see [278] for a 
detailed discussion in the case where the ground state is doubly degenerate. 
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r r 

Figure 20: Lowest energy levels of fc-XORSAT on a c-regular random graph, from exact 
diagonalization. In the inset the region close to the phase transition is magnified. (Left panel) 
Typical instance of the 4-XORSAT model at c = 3, with N = 16. The classical ground state 
degeneracy is 2 4 = 16. (Right panel) USA instance of 3- XORSAT at c = 3, with N = 15. 



is that the first order transition line looks exactly vertical and equal to Tf = 
c/k = 1, suggesting the existence of a hidden duality relating the model at large 
and small T, which was indeed proven in [288, 283]. Note that the spinodal lines 
merge where the first order transition disappears; this point is different from 
the point where the dynamical transition line crosses the first order transition 
line, even if this is not visible in Fig. 18. 

We now discuss the behavior of the gap at the first order transition point 
Tf . There is a definite advantage in this case, namely that a finite fraction 
/usa ~ 0.28 of instances have a single ground state, as discussed above (while 
in previous studies [20, 279, 197] USA instances were exponentially rare). We 
report ED results only on USA instances, in order to unambiguously define the 
gap A(r) between the ground state of H and its first excited state at all values 
of r. The spectrum of a typical USA instance of N = 15 spins is reported 
in Fig. 20. We observe, as expected, that the gap A(r) has a minimum A m j n 
close to the phase transition at Tf (recall that Tf = 1 for c = k at N — > oo). 
In Fig. 21 we show the behavior of the average A m j n as a function of N. Our 
data are clearly consistent with an exponential scaling of the gap, as expected 
based on the discussion of Sec. 4. The probability distribution over instances of 
A m j n has a unique peak close to the average, and its variance is also reported 
in Fig. 21 (dashed bars). This shows that all instances undergo a first order 
transition of the same kind in the thermodynamic limit. These results have 
been confirmed in [283] by computing the minimal gap via QMC (as described 
in Sec. 5.6) at larger sizes. The QMC results confirm the exponential trend. 
Note that in the QMC study the median minimal gap was considered, instead 
of the average. The coincidence of the results confirms that the distribution of 
A m ; n is unimodal and strongly peaked. 

Next, we can compare the cavity results with PIMC. As we already stressed 
several times, in the case c = k = 3, instances have a finite probability of being 
SAT, and otherwise have a minimal energy per spin of order 1/N (see [284, 285]). 
Moreover, a ground state of SAT instances can be found in polynomial time 
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Figure 21: Exact diagonalization data for the minimum gap for fc-XORSAT on c-rcgular 
random graphs. 

(Left panel) k = 4, c = 3, random instances. The classical ground state has an instance- 
dependent degeneracy J\f, so the relevant gap is the one between the ground state and the 
(A/"+l)-th state, see Fig. 20. Full circles represent the average over instances (100 for TV = 8, 12 
and 60 for TV = 16,20) of A m ; n . Full bars represent statistical errors on the average, while 
dashed bars represent the standard deviation over the instances of a single realization of the 
random variable A m i n . Fluctuations are extremely large at small TV, the main contribution 
being due to the fluctuations of N. Dashed line is A min (TV) = 1.244 exp(— 0.065TV) that 
describes well the large TV behavior. 

(Right panel) c = k = 3, USA instances. Full black points represent the average of A m ; n . 
Here ED can be performed up to TV = 24 and fluctuations are reduced: error bars are of the 
order of the symbol size except when explicitly shown (TV = 24). Dashed bars represent the 
standard deviation of a single realization of the random variable A m ; n . Open blue diamonds 
arc QMC data for the median minimal gap from [283]. Here again error bars are of the order 
of symbol size. Dashed line is a fit to A min (TV) = 0.911 exp(-0.081TV). 
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using the Gauss elimination algorithm. This crucial observation allows us to 
find a classical ground state of SAT instances for very large sizes (N = 2049). 
We can therefore run a PIMC starting from the classical ground state at V = 
and slowly increasing T, thereby following the evolution of the classical ground 
states upon introduction of quantum fluctuations. Wc find that the PIMC data 
follow closely the cavity result up to Tf D , see Fig. 18. Then, as expected for 
a first order transition, we find hysteresis around T{ a before the system finally 
jumps to the QP phase. Next, we consider a second PIMC run that is performed 
starting from large r = 2 in the QP phase, and slowly decreasing V. In this case, 
PIMC data follow the cavity ones down to the transition Tf , but then the energy 
remains extensively higher than the ground state energy for any T < T{ . This is 
obviously due to the difficulty in following adiabatically (in the fictitious PIMC 
dynamics) the ground state across an exponentially small gap, as discussed in 
Sec. 5.6.3. This result is an important indication of the difficulty of finding 
the ground state, even in presence of quantum fluctuations, and it is also an 
important proof of the usefulness of the cavity method: in fact, if it were not for 
the Gaussian elimination that allowed us to find the classical ground state and 
run the PIMC starting from it, we would never be able to compute the quantum 
ground state using PIMC. In some models (of which we give an example just 
below) , finding the classical ground states is extremely difficult for any classical 
algorithm. The cavity method allows to compute the ground state energy even 
when Monte Carlo methods fails because of equilibration problems. 

Another demonstration of the difficulty of PIMC to equilibrate (and there- 
fore to find the ground state energy) in this problem is the following. Let us 
consider first the classical limit T = 0. As we discussed in Sec. 3.8, a classical 
Monte Carlo simulation will never be able to equilibrate in polynomial time be- 
low the clustering transition Td(r = 0) [179]. This can be shown by considering 
a classical Monte Carlo simulation that starts in an equilibrium configuration 
at temperature T (which can be generated by the planting technique [169, 196]) 
and computing the spin-spin correlation in Monte Carlo time: 



where the average is done over many realizations of the process. This correlation 
function is reported in Fig. 22. One can see that on approaching the dynamical 
transition the time over which this correlation function goes to zero increases 
as a power law, and it diverges at Td as (T — Td) -7 . Below Td, the correlation 
does not decay to zero anymore, indicating that the dynamics is trapped into a 
state. The same analysis can be repeated at finite T by considering the PIMC 
dynamics. Again, we start by an equilibrium configuration of the paths 8 at 



8 In the quantum case one can still use the planting technique to prepare equilibrium con- 
figurations: the trick consists in doing an initial PIMC run in which both the paths and the 
coupling are changed, i.e. the couplings are treated as dynamical variables. This amounts to 
an "annealed" computation and is equivalent to planting [196]. 




(163) 
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Figure 22: Equilibrium normalized spin-spin correlation functions C(tMCi T)/C(0, T) for the 
Monte Carlo dynamics as functions of tyic f° r several temperatures across the clustering 
transition Td, and averaged over several realizations of the random graph and the initial 
planted configuration. In the inset of each figure, the decorrelation time t(T) such that 
C(t, T) = 1/e (dashed line) is plotted versus T — T d to show the power-law divergence r(T) ~ 
A(T-T d )-~>. 

(Left panel) Classical case T = 0, with standard Metropolis dynamics and N = 60000 spins. 
From left to right, T=0.7, 0.65, 0.6, 0.58, 0.56, 0.55, 0.54, 0.535, 0.53, 0.525, 0.52, 0.515, 0.51, 
0.505. Here T d ~ 0.5098, A ~ 1.82, 7 ~ 3.35. 

(Right panel) Quantum case for T = 0.9, with the PIMC dynamics of [224] and N = 10000. 
From left to right, T=0.45, 0.4, 0.35, 0.33, 0.32, 0.31, 0.30, 0.29, 0.28, 0.27. Here T d ~ 0.303, 
A ~ 3.11, 7 ~ 2.06. The values of T d are consistent with Fig. 18 in both cases. 



temperature T, and we perform a PIMC evolution. We call tyic the PIMC time 
and 

1 r 13 

v*=a dt<Ti{t) (164) 
P Jo 

the imaginary time average of a spin at a given tMC- We define as in the classical 
case 

1 N 

C(t MC ,T) = _^(77 4 (i MC )7J 2 (0)) , (165) 
i=l 

but in this case C(£mc — 0,T) =/= 1 so it is convenient to normalize it by the 
value in £mc = 0. This normalized correlation function is reported in Fig. 22 
and it shows exactly the same behavior as the classical one, with a decorrelation 
time that diverges as (T — Td(r)) _7 when the quantum clustering transition is 
approached. This results shows that PIMC is trapped in a metastable state 
below Td(r) and does not equilibrate. It also shows that Td can be detected via 
a PIMC simulation: this is important because it might allow to estimate Td for 
generic quantum systems for which the cavity method cannot be used, like in 
the classical case [161, 163]. 

6.2.4. UNSAT case: ok 

Here we discuss the UNSAT case c > k, taking as an example c = 4 and 
k = 3. The results for the phase diagram are displayed in Fig. 18. In this case 
the model has a richer phenomenology, very similar to the one of fully connected 
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mean field models [34, 35, 36, 37]. At the RS level, the phenomenology is un- 
changed and a first order is found between the CP and QP phases. However, 
in this case the equilibrium complexity becomes negative below a tempera- 
ture T" c (r) in the dCP phase, signaling that the RS solution becomes incorrect 
(Sec. 5.1.6). The dCP phase then undergoes a thermodynamically second order 
phase transition at T C (T) to a true Spin Glass (SG) phase, where replica sym- 
metry is broken and the Gibbs measure is dominated by a finite number of pure 
states. Therefore, at low enough temperature the first order thermodynamic 
transition happens directly between the SG and QP phases. For this reason, 
the RS computation gives a wrong result for the first order transition line, see 
top panel of Fig. 18. The correct result is obtained by finding the crossing 
between the QP free energy and the SG free energy, and the latter has to be 
computed by optimizing over m the 1RSB free energy and is in general higher 
than the CP free energy as obtained from the RS calculation (Sec. 5.1.6). Still, 
also in this case we conclude on the existence of a first order quantum phase 
transition at T = T c and zero temperature, separating the SG from the QP. The 
transition extends in a line T C (T) ps c/k = 4/3 at low enough temperature, and 
is almost independent of T (at variance with the RS result). 

We tried to repeat the PIMC simulation using the same protocol as in the 
c = k case. However, in this case the problem is typically UNSAT [284, 285]: the 
classical ground states have a finite energy per spin, finding them is extremely 
hard (actually, harder than any NP-complete problem, see Sec. 2.1.2) and the 
quiet planting technique cannot be used. Therefore, in this case we are severely 
limited in the search of the classical ground state, and we can only find it for 
quite small sizes (N < 120) using an exact MAXSAT solver [287]. Still, we could 
repeat the PIMC procedure of increasing T starting from the classical ground 
state, and in this way compute the ground state energy at finite T. A good 
extrapolation in l/N to the thermodynamic limit is possible, and the result 
agrees well with PIQC result, see Fig. 18 and [283]. As in the previous case, we 
find that a PIMC run starting at large T and reducing T fails to find the ground 
state at small T. 

Finally, it is worth to note that the case k = 2 (and any c > k = 2) belongs 
to this class but displays a very different phenomenology. This model has two- 
body interactions and is very close to the Sherrington-Kirkpatrick (SK) model 
(see [152, 153] and Sec. 4), and the transition to the spin glass phase happens 
via a standard second order phase transition instead of a random first order 
transition. In this case, as in the SK model, a second order phase transition 
line in the (T, T) plane separates the CP and SG phases [223] . The transition 
line extends to a quantum critical point at T = 0. In this case there is evidence 
for a polynomially small gap at the critical point [283]. However, as in the 
SK model [208], the spin glass phase seem to be everywhere gapless with an 
exponentially small gap [283]. This should be related to energy- induced level 
crossings inside the spin glass phase, as discussed in Sec. 4.3.2. We refer the 
reader to [283] for a more detailed discussion. 
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6.2.5. Randomization of the transverse fields 

It has been proposed in [41] that small gaps induced by level crossings at 
low values of the transverse field T [39, 40] might be avoided by introducing 
local random fluctuations of the field, i.e. use a different random I\ on each 
spin. Motivated by this proposal, we analyzed whether the first order transition 
can be washed out by a random transverse field. We therefore considered a 
generalization of Eq. (161), where 



€i being a random variable. We choose ej € {1/2,3/2} with equal probability, 
and independently for each spin, in such a way that the average of I\ = Tei 
is equal to T. We solved the model for k = c = 3 using the cavity method 
at the RS level, which we expect to give the exact solution in this case, even 
in presence of random transverse field. The presence of the latter forced us to 
introduce an additional level of population (the model is not factorized anymore, 
see Sec. 5.1.5). We used A4xt = 1000 populations, each containing A/traj = 1000 
trajectories. Again we did not see observable finite population size effects. The 
average free energy as a function of T is qualitatively similar to the one with a 
constant transverse field, showing that the first order phase transition is present 
also with random transverse field. Moreover, because the free energy is self- 
averaging for large N — > oo, a deviation of order 1 in the intensive free energy 
is exponentially rare: this implies that finding a rare sample that does not show 
the transition should be exponentially improbable for large N. We conclude 
therefore that the first order quantum phase transition observed in this model 
is robust against randomization of the transverse field. 

6.2.6. Other approaches 

To conclude this section, we discuss the applicability of the other approaches 
discussed in Sec. 5 to the XORSAT problem. We did not attempt to use vari- 
ational cavity approaches for this problem, because these methods have yet 
to be tested in simpler cases. We discuss here the Operator Quantum Cavity 
(OQC) methods presented in Sec. 5.3. Already at the simplest level £ = 1 (see 
Sec. 5.3.2), the applicability of these methods is severely limited by the need of 
diagonalizing the effective Hamiltonian H c s of Eq. (120). This is because each 
spin in this model interacts with c(k — 1) other spins, therefore H a R is a matrix 
of size 2 c ( fe-1 ) +1 which is quite large already for c = k = 3. For this reason, we 
could only perform a RS variational calculation, which consists in taking opera- 
tor cavity messages at I = 1 of the form (116) with site-independent longitudinal 
and transverse fields, substituting them into the RS free energy, and optimizing 
the latter with respect to the fields. The results are reported in Fig. 23 for 
c = k = 3. Wc found that this approximation gives an extremely good descrip- 
tion of the QP. On the other hand, the description of the CP at large T is not 
excellent: the value of the energy is still an upper bound of the true energy, but 
there is an observable difference at large T, while at the same time the trans- 
verse magnetization is underestimated by the approximation. Moreover, the 
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Figure 23: Comparison between the free energy (right panel) and the transverse magnetization 
(left panel) for XORSAT at k = c = 3, obtained via the path integral quantum cavity method 
(PIQC) and the operator cavity method (OQC). The approximation is almost perfect for the 
QP phase, while it is not very accurate for the CP at large T. 



value of r at spinodal point of the CP is largely overestimated. Despite these 
quantitative discrepancies, we conclude that the variational OQC at I = 1 gives 
a very good qualitative description of the RS phase diagram. Unfortunately, a 
1RSB computation within this approach is out of reach of present computational 
capabilities. 

6.2.7. Discussion 

In this section, we discussed the phase diagram of a typical locked problem: 
the fc-XORSAT model on a c-regular random graph. We have obtained the full 
phase diagram of the quantum model as a function of T and T. The main results 
are: 

1. For k > 2, there is a first order quantum phase transition at T = between 
the low temperature classical phase (which can be either a CP or a SG 
phase) and a QP phase, at a critical value of T = T{ [48]. 

2. The transition is due to a crossing between the low-r classical-like ground 
state(s), and the high-r quantum paramagnetic state. It is of very differ- 
ent nature from the level crossing at small T between different spin glass 
ground states discussed in [39, 40, 41]. 

3. The first order transition is observed for almost all instances, even for very 
small N. In general, finite size effects are extremely small in this model, 
and they are mainly due to the fluctuations in the number of classical 
ground states. 
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4. The first order transition is generically associated to an exponentially van- 
ishing gap of H [38, 48], hence, in this model, the QAA requires a run 
time scaling exponentially with the system size to find the ground state. 

5. The case k = 2 is special: here the model has a second order transition 
with a polynomially small gap, but level crossings at small T are present 
and induce an exponentially small gap, hence also in this case the QAA 
is not efficient [283]. 

The main missing ingredient in locked models with respect to the general picture 
outlined in Sec. 3.6 is the internal entropy of the clusters. In these models 
clusters are isolated configurations, and they are very far away from each other. 
If the local environment around each solution is the same, as in the XORSAT 
problem, the degeneracy is not lifted at any order in perturbation theory. Then, 
the degeneracy is only lifted by (non-pcrturbative) quantum tunneling leading 
to an exponentially small splitting. Level crossings at small T are therefore 
very difficult to observe in these models, and in the thermodynamic limit the 
ground states remain closely degenerate on increasing T up to the first order 
phase transition. 

6.3. The coloring problem 

In order to discuss the effect of the existence of exponentially many clus- 
ters (pure states) with a non-trivial distribution of internal energy and entropy 
densities, we consider here the quantum version of the simplest classical model 
that shows this effect: the coloring problem on random regular graphs. The 
phenomenology of this problem is rather different from the one of XORSAT: 
contrary to the latter, the coloring problem is not a locked model. Hence, in the 
classical limit the number of solutions is exponentially large in N in the SAT 
phase. As it has been discussed in Sec. 3.7, the structure of the solution space 
of this model is very similar to the one of the Random Subcubes Model (RSM) 
discussed in Sec. 3.6. We expect that adding quantum fluctuations should lead 
to an extremely complex spectrum, as it has been discussed in Sec. 4.4 for the 
Quantum RSM (QRSM). In particular, we want to show that due to entropic 
effects, different clusters evolve very differently under the action of the quantum 
term, leading to level crossings. 

6.3.1. Definition of the model 

The model is defined as follows. We consider N Potts spins. Each of these 
is a classical variable <7j S X = {1, • • • , q}. To define the quantum model we 
introduce for each spin the Hilbert space spanned by \o~i) and the operators 




(167) 
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where 5 aiCT i is the Kronecker delta. We define the Hamiltonian as 

where the first sum runs on pairs of variables that are connected by a link 
of the random regular graph with fixed connectivity c, on which the model is 
defined. In the classical case T = it is easy to check that the model reduces 
to the classical coloring problem (or equivalently to the antiferromagnetic Potts 
model) introduced at the beginning of Sec. 2.1. In the following we will call 
m x = (Tj) the "transverse magnetization" by analogy with the Ising spin case. 

We want to stress that the analysis of this model turns out to be extremely 
difficult for a combination of technical reasons. First of all, the interesting 
regime is q > 4 (the case q = 3 being special, because it displays a continuous 
transition at the classical level) and quite large connectivities c [173]. Hence, 
exact diagonalization is impossible, because the size of the Hilbert space grows 
as q N , much faster than for Ising spins, and moreover for large c one needs to 
consider quite large systems to avoid finite size effects. One therefore cannot 
have any direct information about the spectrum and needs to infer it from other 
techniques. OQC methods are impossible to use, again because the size of the 
Hilbert space grows too quickly with the number of spins. 

We mainly used the PIQC method, which however is more difficult to apply 
to this model than to the XORSAT model. We found that large populations 
have to be considered to avoid finite population size effects: for the 1RSB com- 
putations reported below, we typically used A/jut ~ 4000 populations made of 
A4raj ~ 4000 trajectories, as in the XORSAT case (the coloring problem is also 
factorized). However, for large q and c the PIQC solution algorithm is slow (the 
running time grows roughly as q 2 c, see the discussion in Sec. 5). Moreover, in 
many cases we had to perform finite population size scalings to obtain reliable 
results. We also used PIMC, which requires similar computational effort than 
in other models, but suffers as usual from equilibration problems. 

Keeping these difficulties in mind, we now proceed to discuss briefly the 
structure of this model. We will keep the discussion short, and we will focus 
on the computationally simplest case that displays the phenomenology we are 
interested in, namely q = 4 and c = 9. A much more complete and detailed 
discussion of this problem can be found in [49] . 

6.3.2. Results 

For q = 4 and c = 9 the model is classically satisfiable (the ground state 
energy is zero and graphs are typically colorable). At T = the number of 
ground states (solutions) is exponentially large, and they are arranged in an 
exponential number of clusters [173]. This corresponds to the region Cd < c < c c 
according to the discussion of Sec. 3. In other words, T c = while Td = 0.153(5) 
for the classical model at T = 0, which is described by the 1RSB cavity equations 
at Parisi parameter m = 1, below Td and down to T = 0. 
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Let us briefly recall the expected behavior for small T > and low tem- 
perature T, based on the analysis of the QRSM (Sec. 4.4). At T = 0, clusters 
have internal entropies that are distributed according to a large deviation func- 
tion (the complexity), Af(s) = exp[iV£(s)]. Typical configurations are found 
in clusters that have a value of the entropy s* such that E'(s*) = —1, and 
the corresponding complexity S(s*) is strictly positive. However, many other 
clusters with larger and smaller entropies exist. When r > 0, each cluster A of 
degenerate classical states transforms continuously into a set of quantum states, 
the lowest of which (the "ground state of cluster A" \GS(A))) has an energy 
per spin 

e(r) =e cl (A,T)-Tm x (A,T) , 
e cl (A, T) = (GS(A)\Hp\GS(A))/N , (lg9) 
m x (A,T) = (GS(A)\Y,Ti\GS(A))/N . 

i 

The crucial observation is that, as in the QRSM, we expect m x (A,T) to be 
finite when r — > 0, limr-s-o m-x(A, T) — m x (A), and we expect m x (A) to be 
positively correlated with the classical entropy of the cluster. At the same time, 
e c i(A,r) oc T 2 at small T. 

Therefore, the energy of the ground state of a cluster is linear at small T, 
e(r) ~ — Tm x (A). Largest clusters yield the greatest decrease in energy when 
quantum fluctuations are switched on, and they dominate at zero temperature 
as soon as L > 0. Because these are the states with maximal entropy, they 
correspond to £(s max ) = 0. Hence as soon as T > 0, the zero temperature 
complexity drops abruptly to zero. In other words, we expect the system to 
condense into the largest cluster under an infinitesimal amount of quantum 
fluctuations. This in particular implies that a non-zero T C (T) should emerge, 
and that T c cx T for small T, as in the QRSM. 

This scenario is indeed confirmed by the PIQC results, sec figure 24. We 
find that starting from finite temperature and increasing L, the equilibrium 
complexity at Parisi parameter m = 1 decreases from its (finite) classical value, 
until it vanishes at some T C (T). Inverting this function we find T C (T), which 
is indeed a linearly increasing function of T at small T. This result is a strong 
indirect confirmation that the QRSM describes correctly the physics of complex 
models like the quantum coloring. Next, we examine the evolution with F of 
7d, which is found by performing scans at constant T and m = 1, starting at 
low T and increasing T until the non-trivial 1RSB solution is lost. We find that 
also Td increases slightly from its classical value for small F > 0. The behavior 
of and T c shows that, contrary to naive intuition, quantum fluctuations 
promote glass formation, rather than make the glass unstable. This has been 
recently found in a series of different models of glasses by mean of different 
methods [47, 198, 281, 289]. 

At larger F, the two lines 7d(r) and T C (T) approach each other and at some 
value of r they cross a third line Ti(r) (see Fig. 24), that corresponds to a 
linear instability of the RS solution of the cavity equations towards 1RSB (i.e. 
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Figure 24: Phase diagram of the quantum coloring problem for q = 4 and c = 9. The black 
solid line with circles represents the linear instability of the RS solution, T;(r) (a linear fit 
at small T is plotted with a black dashed line); the blue line with squares is the clustering 
transition T(j(r) separating the P and dP phases; the green solid line with triangles is the 
condensation transition T c (r) separating the dP and SG phases (dashed green line is a linear 
fit of the small T behavior). The classical problem at T = is in the SAT phase, with an 
exponential number of clusters. It has therefore a finite T<j ~ 0.f5 and T c = 0. For small 
r > 0, we observe a linear increase of T C (T) ~ T, as in the RSM, confirming that the physics 
of the two models is very similar. At larger T the transition becomes a continuous transition 
Ti(r) that extends to a third order quantum critical point at Ti ~ 0.47. 

a continuous RSB transition). Because this instability is a property of the RS 
solution, it is independent of m. Hence it has been detected by solving the 
1RSB cavity equations at m = 0, initializing the population at T = and T > 
in the RS solution (all external populations are identical), increasing Y and 
finding the point where a 1RSB solution appears continuously. When the lines 
cross, the transition ceases to be a discontinuous 1RSB transition (i.e. a random 
first order transition) and becomes instead a continuous RSB transition. We 
find that at larger Y the transition remains continuous until T[(Y) goes to zero 
around Yi ~ 0.45, leading to a second order quantum phase transition (third- 
order in the thermodynamic sense), at variance with the XORSAT case where 
the transition is first order. Correspondingly, in this case the paramagnetic 
phase is unique and we refer to it as P in Fig. 24. 

In Fig. 25 we show the energy and transverse magnetization as a function 
of r at fixed temperature T = 0.06, as obtained from PIQC and PIMC. For 
this temperature, the system undergoes first a condensation transition at Y c ~ 
0.0427 (from dP to SG), then a second order transition at Yi ~ 0.45 (from SG 
to P). The former transition is of second order, while the latter is of third order. 
They both lead to weak singularities in the derivatives of e and m x that are not 
visible in the figure. The PIQC computations arc 1RSB at m = 1 for Y < Y c , 
while m — m* < 1 for T c < Y < Y[. Above T;, the RS computation is correct. 

PIMC simulations have been performed in two ways. In the first run, we 
prepared a typical graph together with one of its typical solutions via the "quiet 
planting" technique [196]; we initialized the PIMC at Y = in the solution, 
hence at zero energy, and we then slowly increased Y. In the second run, we ini- 
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Figure 25: Energy (left panel) and transverse magnetization (right panel) as a function of T 
at fixed temperature T = 0.06 for q = 4 and c = 9. PIQC computations are reported as black 
solid lines. PIMC simulations (for a single system of size TV = 10000) are reported as a dashed 
black line with squares for increasing T and as a dashed red line with circles for decreasing T. 
(Inset of left panel) Zoom on the region of low T where hysteresis is observed in the PIMC. 
The blue arrow indicates the value of classical energy (0.0019) that corresponds to a classical 
annealing starting from Td. 

tialized the PIMC at T = 2, we equilibrated in the paramagnetic phase, and we 
decreased T down to T = 0. On the scales of the figure, no difference between the 
two PIMC runs is observed, however a closer look (inset of left panel in Fig. 25) 
reveals that decreasing T one obtains a positive residual energy at T = 0. The 
latter is found to be larger than the residual energy after an infinitely slow 
thermal annealing, which is obtained by preparing a quietly planted configura- 
tion at Td and performing a slow classical annealing down to T = [196, 187]. 
Although a direct comparison is not possible because the PIMC annealing has 
not been extrapolated in the infinitely slow limit, this result suggests that an 
annealing of an imaginary time PIMC simulation (Sec. 5.6.3) is not more effi- 
cient than a thermal annealing for this model, and we believe that this is once 
again related to entropic level crossings inside the SG phase, as in the QRSM 
model. A more detailed investigation of this point is one of the most important 
directions for future research. 

Finally, we want to confirm numerically the assumptions in Eq. (169). Un- 
fortunately, neither PIQC nor PIMC can access directly the zero temperature 
properties of a cluster. To solve the problem, we performed several PIMC runs at 
different low temperatures, starting from the same quietly planted zero-energy 
configuration at T = and increasing T. In this way we assume that the PIMC 
follows the evolution of a single cluster, selected by the planted configuration, 
in T and T. Extrapolating the results to T — > gives the ground state proper- 
ties of the cluster. In the QRSM, the transverse magnetization of a cluster is 
m x (A) = s(A) tanh(/?r), hence it vanishes linearly in T at any finite temper- 
ature, but with a diverging slope that indicates a finite m x in the limit where 
T — > first, and then T —> 0. Note that also in the coloring problem one can 
show in perturbation theory that whenever floppy spins (that can be flipped 
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Figure 26: PIMC results for a single instance with q = 4, c = 9 and TV = 1000 spins, starting 
from a planted state at T = and increasing T at fixed T. The thick black line is the 
extrapolation at T = 0, that shows that limp_>o m x (A, T) = mj > 0. 

without energy change) are present in a cluster, the slope of m x computed by 
perturbation theory for this cluster diverges when T — > [49]. The results of 
the PIMC runs, reported in Fig. 26, are consistent with this expectation and the 
extrapolation clearly shows a finite m x at T = 0. In order to complete 

this investigation, one should show that larger clusters have larger m x ; however 
generating configurations belonging to clusters of (appreciably) different size is 
not possible using the quiet planting, leaving a direct investigation of this point 
as an open problem. 

6.3.3. Discussion 

Compared to the XORSAT model, the coloring problem shows a much more 
complex phase diagram, due to the non-trivial distribution of the energy and 
entropy densities of its pure states. We studied the problem for q = 4 and c = 9, 
that classically is in the clustered phase with an exponential number of clusters. 
By means of combined PIQC and PIMC computations, we showed that: 

• The zero temperature transition from the spin glass to the paramagnetic 
phase is of third order, unlike in the XORSAT problem. Although we 
suspect that the transition becomes first order for large enough q and c, 
we cannot investigate the problem because of computational limitations. 

• The quantum spin glass phase has a complex structure similar to the one 
of the QRSM. Due to the finite entropy of clusters, the states are ex- 
tremely delocalized even for T — > 0, leading at T = to a finite transverse 
magnetization m x for r — > and a linearly decreasing ground state energy 
e(r) ~ -Tm° x . 
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• We cannot give direct evidence that m° grows with the classical entropy of 
the cluster. However, this is indirectly confirmed by the fact that T c rises 
linearly from zero for small T, as in the QRSM. Because of this entropic 
phenomenon, quantum fluctuations promote the formation of the glass at 
low enough temperature and T, as found in [47, 198, 281, 289]. Hence, we 
expect that level crossings are induced in this model by an energy-entropy 
effect, exactly as in the QRSM. We believe that these crossings will lead 
to an everywhere exponentially small gap in the spin glass phase. 

We expect that similar results about level crossings hold for any other random 
optimization problem characterized by clusters of finite entropy However, the 
order of the zero temperature transition depends on the problem under investi- 
gation. 

7. Conclusions 

Let us now summarize the main messages we wanted to convey in this review, 
and draw some perspectives for future work. As the Quantum Adiabatic Algo- 
rithm (QAA) is a general purpose optimization procedure, random constraint 
satisfaction problems provide natural benchmarks for its efficiency The tools 
of the statistical mechanics of disordered systems that have been developed for 
their study, first in the classical case and more recently in a quantum context, 
have several main outcomes. 

Natural benchmarks for the QAA are random instances with a Unique Sat- 
isfying Assignment (USA), because for these instances the minimal gap in un- 
ambiguously defined. The generation process of USA instances is enlightened 
by the detailed knowledge acquired on classical random constraint satisfaction 
problems; most of these ensembles display a satisfiability transition that is 
strongly discontinuous in terms of the number of their solutions, which jumps 
from being exponentially large to zero when the control parameter is tuned 
across the transition. In consequence USA have an exponentially small prob- 
ability of existence even at the transition, so that alternative ensembles where 
USA exist with finite probability should instead be used (Sec. 3.9). 

The classical phase diagrams of random constraint satisfaction problems ex- 
hibits a rich phenomenology with several structural phase transitions, that are 
not all present in every model (Sec. 3). When quantum fluctuations arc added 
zero temperature phase transitions as a function of the intensity of the fluctua- 
tions gcnerically appear. In some models they are of the first order type, thus 
accompanied by an exponentially small gap that implies an exponentially large 
running time for the QAA. Other models exhibit higher order phase transitions, 
so that the gap at the transition can be only polynomially small, yet in general 
their weakly quantum phase is gapless because of the complex spin glass struc- 
ture that induces a continuum of avoided level crossings through the competition 
between the energy and the entropy of the classical pure states (Sec. 6). The 
path integral quantum cavity method provides an analytic framework in which 
these transitions can be quantitatively computed in the thermodynamic limit. 
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Even if it contains an heavy numerical part, corresponding to the resolution of 
the cavity equations, it is of a different nature compared to quantum Monte 
Carlo methods that are plagued by equilibration problems when glassy features 
of random constraint satisfaction problems come into play. It is also worth to 
note that the methods discussed in Sec. 5 are useful also for condensed matter 
applications to bosonic or fermionic systems (see e.g. [232, 290, 229, 198]). 

We would also like to emphasize that, as discussed for instance in Sec. 3.8, 
there always exists an annealing (be it classical or quantum) path that allows 
to find a solution in polynomial time, obtained by adding an external field in 
the direction of one solution of the problem. The relevant question, in our 
opinion, is not about proving the existence of good annealing paths, but to 
determine whether an efficient annealing path can be constructed without as- 
suming a detailed knowledge of the ground states of the Hamiltonian to be 
minimized. In this review we mainly considered annealing in a (possibly ran- 
dom) transverse field, but we believe that our results can be applied to more 
general Hamiltonians, at least those which do not have a sign problem, e.g. 
Bosonic ones [232, 290, 198]. 

Even though the above statement sounds rather negative for the usefulness 
of the quantum adiabatic algorithm, we want to emphasize that it might be 
more promising to consider it as an approximation algorithm and not as an 
exact one. Classical computational complexity theory has indeed established 
several hardness of approximation results and it would be interesting to study 
whether a quantum annealing performed on a timescale that does not respect the 
adiabaticity condition could lead to energies that, even if higher than the ground 
state one, would still be lower than the best approximation ratio achievable 
by classical approximation algorithms. This question is obviously much more 
difficult to tackle mathematically, as there are no general results as the adiabatic 
theorem to bound the residual energy after a non-adiabatic evolution. 

Finally, an important aspect that has not been touched upon at all in this 
review is the description of actual quantum computers, that should necessarily 
take into account the unavoidable coupling to the environment. Quantum adia- 
batic algorithm is claimed to be robust with respect to a thermal coupling with 
the environment [291, 292, 293]; we believe it would be worth investigating the 
modification of the picture discussed in this review induced by an experimentally 
relevant modelization of the coupling with the environment. 
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